Intercerebral current source estimation method intercerebral current source estimation program recording medium containing the intercerebral current source estimation program and intercerebral current source estimation apparatus

ABSTRACT

Based on the fact that by causing an appropriate current flow on a virtual curved surface between a current source and an observation surface, an electromagnetic field generated by the current source can be recovered, and that as the curved surface comes closer to the true current source, current distribution on the curved surface becomes smaller, Bayesian estimation of the current source that recovers the observed data is performed. In this estimation, the fact that the model posterior probability attains the maximum when the curved surface includes the current source is utilized, that is, the model posterior probability is considered, to estimate the position of the current source including the depth direction. When observation data obtained by other observation means is available simultaneously, such information is incorporated in hierarchical prior distribution for Bayesian estimation, to enable estimation with higher precision.

TECHNICAL FIELD

The present invention relates to a method of estimating position ordistribution of a brain current source or sources generating magneticfields or electric fields outside a scalp in correspondence to brainactivities, as well as to a brain current source estimating program, arecording medium recording the brain current source estimating programand to a brain current source estimating apparatus.

BACKGROUND ART

Thanks to remarkable development in the field of biomedical measurementtechnique, precision in measuring weak electric field (brain wave) orweak magnetic field (brain magnetic wave) generated from the brain, ofwhich measurement has been difficult and error-prone, has been improvedyear by year.

Specifically, receiving external stimuli, neural cells in the braingenerate a current. The current results in the afore-mentioned weakelectric field or weak magnetic field. Here, “brain wave” refers to theelectric field generated from the brain by the current from the neuralcells. An “electroencephalogram: EEG” refers to a method of measuringthe brain wave.

The “brain magnetic wave” refers to a magnetic field generated from thebrain by the current from the neural cells. A “magnetoencephalography:MEG” refers to a method of measuring the brain magnetic wave. A crucialadvantage of the magnetoencephalography is that, as the magnetic fieldis almost free of any influence of volume conductor, it is expected thatrelatively accurate three-dimensional estimation of the position of abrain current source can be attained by measuring magnetism from outsideone's scalp.

In the analysis of the brain magnetic wave, an active portion of thebrain is estimated in a non-invasive manner, by measuring the generatedmagnetic field from outside the brain.

The magnetic field, however, is so weak that it is very susceptible tothe influence of external magnetic field such as terrestrial magnetism.Therefore, the weak magnetic field is measured by a SuperconductingQUantum Interference Device (SQUID), which is a measuring deviceutilizing superconductivity, within a shield that shuts out any externalmagnetic field.

It is noted, however, that in the field of studying algorithms for“estimating positions of brain current source,” decisive method isnon-existent at present, though various and many variations of initialmodels have been tried.

By way of example, “dipole estimation method” as one algorithm for“estimating positions of brain current sources,” is disclosed inReference 1: J. C. Mosher, P. S. Lewis and R. M. Leahy, IEEE Trans.Biomed. Engng. <1992> vol. 39, pp.541-557. In the “dipole estimationmethod,” however, the position of a dipole is estimated from observedmagnetic field, assuming that the current source in the brain can berepresented by one or a number of current dipoles, and this method isdisadvantageous in that it is difficult to determine the number ofdipoles.

As another algorithm, “spatial filtering method” is disclosed inReference 2: K. Toyama, K. Yoshikawa, Y. Yoshida, Y. Kondo, S. Tomita,Y. Takanashi, Y. Ejima and S. Yoshizawa, Neuroscience, 1999, 91 (2), pp.405-415. In the “spatial filtering method,” location of a brain currentsource is restricted in consideration of physiological findings, anddistribution of dipoles is estimated. This method is disadvantageous inthat accurate estimation of the depth of the current source is notpossible.

DISCLOSURE OF THE INVENTION

An object of the present invention is to provide a brain current sourceestimating method that enables estimation of a position, depth directioninclusive, of a brain current source.

Another object of the present invention is to provide a brain currentsource estimating method that enables estimation of positions of aplurality of brain current sources from observed data, even when thereare a plurality of brain current sources.

A still further object of the present invention is to provide a braincurrent source estimating method enabling estimation with higheraccuracy, when observation data obtained by a method of observationindependent of the observation of magnetic field for estimating thebrain current source are available, by combining data obtained by theplurality of observation methods.

A still further object of the present invention is to provide a braincurrent source estimating program that enables estimation of a position,depth direction inclusive, of a brain current source and a recordingmedium having the program recorded thereon.

A still further object of the present invention is to provide a braincurrent source estimating program that enables estimation of positionsof a plurality of brain current sources from observed data, even whenthere are a plurality of brain current sources, and a recording mediumhaving the program recorded thereon.

A still further object of the present invention is to provide a braincurrent source estimating program enabling estimation with higheraccuracy, when observation data obtained by a method of observationindependent of the observation of magnetic field for estimating thebrain current source are available, by combining data obtained by theplurality of observation methods, and a recording medium having theprogram recorded thereon.

A still further object of the present invention is to provide a braincurrent source estimating apparatus that enables estimation of aposition, depth direction inclusive, of a brain current source.

A still further object of the present invention is to provide a braincurrent source estimating apparatus that enables estimation of positionsof a plurality of brain current sources from observed data, even whenthere are a plurality of brain current sources.

A still further object of the present invention is to provide a braincurrent source estimating apparatus enabling estimation with higheraccuracy, when observation data obtained by a method of observationindependent of the observation of magnetic field for estimating thebrain current source are available, by combining data obtained by theplurality of observation methods.

In order to attain the above described objects, the present inventionprovides a brain current source estimating method for estimating, basedon an electromagnetic field observed outside a scalp, a position of acurrent source as a source of the electromagnetic wave existing in thebrain, including the steps of: setting, in the brain, a plurality ofvirtual curved surfaces having depths from brain surface different fromeach other and shapes not intersecting with each other and settinglattice points on each of the virtual curved surfaces; estimating, oneach of the virtual curved surfaces, a current distribution forrecovering the observed electromagnetic field; and based on an expansionof the current distribution estimated on the virtual curved surfaces anda difference between the electromagnetic field recovered based on thecurrent distribution and the observed electromagnetic field, identifyinga virtual curved surface at which the expansion and the differenceattain relative minimums among the plurality of virtual curved surfacesas a true curved surface on which the current source exists.

Preferably, the step of estimating the current distribution includes thestep of determining posterior probability by Bayesian estimation methodfrom prior distribution and observation data of the electromagneticfield; and the step of identifying as a true curved surface on which thecurrent source exists includes the step of identifying a virtual curvedsurface of which corresponding posterior probability attains thehighest, among the virtual curved surfaces.

Preferably, the step of estimating a current distribution includes thestep of identifying a first virtual curved surface closest to the brainsurface and having posterior probability attaining a relative maximum,among the plurality of virtual surfaces, while successively moving froma virtual curved surface on the side of the brain surface to a deeperside; and the step of identifying a curved surface as a true curvedsurface on which the current source exists includes the steps ofidentifying a localized current distribution corresponding to a point ofrelative maximum of the current distribution, on the first virtualcurved surface, separating a plurality of local surfaces each includingthe localized current distribution, and fixing, among the plurality oflocal surfaces, local surfaces other than a local surface as an objectof identification, moving the local surface as an object ofidentification in the depth direction, and identifying positions wherethe posterior probability attains the relative maximum, successivelyfrom the side closer to the brain surface.

Preferably, in the step of estimating a current distribution, thecurrent distribution is estimated with a first spatial resolution; andthe method further includes the step of re-estimating the currentdistribution with a second spatial resolution higher than the firstresolution and resolution of the plurality of virtual curved surfaces inthe depth direction being improved.

Preferably, the step of estimating a current distribution includes thestep of setting, when the current distribution is estimated inaccordance with Bayesian estimation, a hierarchical prior distributionin the Bayesian estimation using observation data obtained by otherobservation method independent of the observation of electromagneticfield for the estimation of the current source. By way of example, whenit is known from observation data obtained by a different method ofobservation that the location of the brain current source is within arestricted area, search in the depth direction may be omitted and thecurrent distribution in the restricted area can be obtained by Bayesianestimation.

According to another aspect, the present invention provides a programfor a computer for estimating, based on an electromagnetic fieldobserved outside a scalp, a position of a current source as a source ofthe electromagnetic wave existing in the brain, to have the computerexecute the steps of: setting, in the brain, a plurality of virtualcurved surfaces having depths from brain surface different from eachother and shapes not intersecting with each other and setting latticepoints on each of the virtual curved surfaces; estimating, on each ofthe virtual curved surfaces, a current distribution for recovering theobserved electromagnetic field; and based on an expansion of the currentdistribution estimated on the virtual curved surfaces and a differencebetween the electromagnetic field recovered based on the currentdistribution and the observed electromagnetic field, identifying avirtual curved surface at which the expansion and the difference attainrelative minimums among the plurality of virtual curved surfaces as atrue curved surface on which the current source exists.

Preferably, the step of estimating the current distribution includes thestep of determining posterior probability by Bayesian estimation methodfrom prior distribution and observation data of the electromagneticfield; and the step of identifying as a true curved surface on which thecurrent source exists includes the step of identifying a virtual curvedsurface of which corresponding posterior probability attains thehighest, among the virtual curved surfaces.

Preferably, the step of estimating a current distribution includes thestep of identifying a first virtual curved surface closest to the brainsurface and having posterior probability attaining a relative maximum,among the plurality of virtual surfaces, while successively moving froma virtual curved surface on the side of the brain surface to a deeperside; and the step of identifying a curved surface as a true curvedsurface on which the current source exists includes the steps ofidentifying a localized current distribution corresponding to a point ofrelative maximum of the current distribution, on the first virtualcurved surface, separating a plurality of local surfaces each includingthe localized current distribution, and fixing, among the plurality oflocal surfaces, local surfaces other than a local surface as an objectof identification, moving the local surface as an object ofidentification in the depth direction, and identifying positions wherethe posterior probability attains the relative maximum, successivelyfrom the side closer to the brain surface.

Preferably, in the step of estimating a current distribution, thecurrent distribution is estimated with a first spatial resolution; andthe method further includes the step of re-estimating the currentdistribution with a second spatial resolution higher than the firstresolution and resolution of the plurality of virtual curved surfaces inthe depth direction being improved.

Preferably, the step of estimating a current distribution includes thestep of setting, when the current distribution is estimated inaccordance with Bayesian estimation, a hierarchical prior distributionin the Bayesian estimation using observation data obtained by otherobservation method independent of the observation of electromagneticfield for the estimation of the current source. By way of example, whenit is known from observation data obtained by a different method ofobservation that the location of the brain current source is within arestricted area, search in the depth direction may be omitted and thecurrent distribution in the restricted area can be obtained by Bayesianestimation.

According to a still further aspect, the present invention provides abrain current source estimating apparatus for estimating, based on anelectromagnetic field observed outside a scalp, a position of a currentsource as a source of the electromagnetic wave existing in the brain,including: virtual curved surface setting means for setting, in thebrain, a plurality of virtual curved surfaces having depths from brainsurface different from each other and shapes not intersecting with eachother and setting lattice points on each of the virtual curved surfaces;current distribution estimating means for estimating, on each of thevirtual curved surfaces, a current distribution for recovering theobserved electromagnetic field; and current source identifying means foridentifying, based on an expansion of the current distribution estimatedon the virtual curved surfaces and a difference between theelectromagnetic field recovered based on the current distribution andthe observed electromagnetic field, a virtual curved surface at whichthe expansion and the difference attain relative minimums among theplurality of virtual curved surfaces as a true curved surface on whichthe current source exists.

Preferably, the current distribution estimating means includes posteriorprobability determining means for determining posterior probability byBayesian estimation method from prior distribution and observation dataof the electromagnetic field; and the current source identifying meansincludes virtual curved surface identifying means for identifying avirtual curved surface of which corresponding posterior probabilityattains the highest, among the virtual curved surfaces.

Preferably, the current distribution estimating means includesshallowest virtual curved surface identifying means for identifying afirst virtual curved surface closest to the brain surface and havingposterior probability attaining a relative maximum, among the pluralityof virtual surfaces, while successively moving from a virtual curvedsurface on the side of the brain surface to a deeper side; and thecurrent source identifying means includes localized current distributionidentifying means for identifying a localized current distributioncorresponding to a point of relative maximum of the currentdistribution, on the first virtual curved surface, local surfaceextracting means for separating a plurality of local surfaces eachincluding the localized current distribution, and local surface positionidentifying means for fixing, among the plurality of local surfaces,local surfaces other than a local surface as an object ofidentification, moving the local surface as an object of identificationin the depth direction, and identifying positions where the posteriorprobability attains the relative maximum, successively from the sidecloser to the brain surface.

Preferably, the current distribution estimating means estimates thecurrent distribution with a first spatial resolution and thereafterre-estimates the current distribution with a second spatial resolutionhigher than the first resolution and resolution of the plurality ofvirtual curved surfaces in the depth direction being improved.

Preferably, the current distribution estimating means includes means forsetting, when the current distribution is estimated in accordance withBayesian estimation, a hierarchical prior distribution in the Bayesianestimation using observation data obtained by other observation methodindependent of the observation of electromagnetic field for theestimation of the current source. By way of example, when it is knownfrom observation data obtained by a different method of observation thatthe location of the brain current source is within a restricted area,search in the depth direction may be omitted and the currentdistribution in the restricted area can be obtained by Bayesianestimation.

According to the brain current source estimating method of the presentinvention, it is possible to estimate the position of a brain currentsource including the depth direction, from observation data of MEG orEEG. Such estimation in the depth direction is applicable even whenthere are a plurality of current sources. Further, the method isapplicable where the current sources exist locally as in the case ofcurrent dipoles and applicable to a current source that has a certainexpansion. Further, it is possible to estimate how the current sourceexpands.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic illustration of an exemplary configuration of anMEG system.

FIG. 2 is a schematic illustration representing a magnetic fieldgenerated by a current source, observed on an appropriate curvedsurface.

FIG. 3 is a schematic illustration representing a relation between acurrent source in the brain and a plurality of virtual curved surfaces.

FIG. 4 is a flow chart representing an overall procedure of the braincurrent source estimating method in accordance with the presentinvention.

FIG. 5 is a flow chart representing a process of initial estimation ofthe current source using an initial resolution.

FIG. 6 is a flow chart representing a process for identifying a currentsource closest to the brain surface.

FIG. 7 is a first flow chart representing a process for identifyingdepth of a current source corresponding to each local surface.

FIG. 8 is a second flow chart representing a process for identifyingdepth of a current source corresponding to each local surface.

FIG. 9 is a first flow chart representing a process for re-estimating aposition of a current source with higher spatial resolution.

FIG. 10 is a second flow chart representing a process for re-estimatinga position of a current source with higher spatial resolution.

FIG. 11 is a top view of a magnetic field distribution observed on asurface of a hemisphere, assuming that the human brain is a hemisphere.

FIGS. 12(a), 12(b) and 12(c) represent results of initial estimation ofcurrent sources using initial resolution, where the radius is R=5.0,R=6.0 and R=7.0, respectively.

FIGS. 13(d), 13(e) and 13(f) represent results of initial estimation ofcurrent sources using initial resolution, where the radius is R=7.0,R=8.0 and R=9.0, respectively, that is, closer to the surface of thebrain.

FIG. 14 represents magnitude of free energy on each virtual hemisphere,obtained when current distribution that attains highest free energy iscalculated.

FIGS. 15(a), 15(b) and 15(c) represent processes for identifying thecurrent source executed further on a local surface including a point ofrelative maximum, where the radius is R=5.0, R=6.0 and R=7.0,respectively.

FIGS. 16(d), 16(e) and 16(f) represent processes for identifying thecurrent source executed further on a local surface including a point ofrelative maximum, where the radius is R=7.5, R=8.0 and R=9.0,respectively.

FIG. 17 represents the free energy calculated with the depth of localsurface varied.

FIG. 18 is a top view of a magnetic field distribution observed on asurface of the human brain assumed as a hemisphere, when there are twocurrent sources in the brain.

FIGS. 19(a), 19(b) and 19(c) represent results of initial estimation ofcurrent sources using initial resolution, where the radius is R=5.0,R=6.0 and R=7.0, respectively.

FIGS. 20(d), 20(e) and 20(f) represent results of initial estimation ofcurrent sources using initial resolution, where the radius is R=7.5,R=8.0 and R=9.0, respectively.

FIG. 21 represents radius dependency of free energy, when currentdistribution that attains highest free energy is calculated for eachvirtual hemispherical surface.

FIGS. 22(a), 22(b) and 22(c) represent current distribution on localsurfaces where the radius is R=5.0, R=6.0 and R=7.0, respectively,illustrating the process for calculating the depth of a first localsurface attaining maximum posterior probability, to identify a currentsource closest to the brain surface.

FIGS. 23(d), 23(e) and 23(f) represent current distribution on localsurfaces where the radius is R=7.5, R=8.0 and R=9.0, respectively,illustrating the process for calculating the depth of a first localsurface attaining maximum posterior probability, to identify a currentsource closest to the brain surface.

FIG. 24 represents local surface position (radius) dependency of freeenergy when a virtual local surface is moved.

FIGS. 25(a), 25(b) and 25(c) represent current distribution on localsurfaces where the radius is R=5.0, R=5.5 and R=6.0, when one localsurface is fixed on a specific surface and a local surface correspondingto the other current source is moved.

FIGS. 26(d), 26(e) and 26(f) represent current distribution on localsurfaces where the radius is R=6.5, R=7.0 and R=7.5, when one localsurface is fixed on a specific surface and a local surface correspondingto the other current source is moved.

FIG. 27 represents local surface position (radius) dependency of freeenergy when a virtual local surface is moved.

FIG. 28 is a first graph representing a result of simulation consideringlocalized condition only.

FIG. 29 is a second graph representing a result of simulationconsidering localized condition only.

FIG. 30 is a first graph representing a result of simulation consideringboth continuous condition and localized condition.

FIG. 31 is a second graph representing a result of simulationconsidering both continuous condition and localized condition.

FIG. 32 is a first graph representing a result of simulation consideringboth continuous condition and localized condition, andusing-hierarchical prior distribution that facilitates estimation of acurrent in a region where presence of a current source is highly likely.

FIG. 33 is a second graph representing a result of simulationconsidering both continuous condition and localized condition, and usinghierarchical prior distribution that facilitates estimation of a currentin a region where presence of a current source is highly likely.

BEST MODES FOR CARRYING OUT THE INVENTION

Embodiments of the present invention will be described with reference tothe figures.

As already described, magnetoencephaography: MEG andelectroencephalogram: EEG have been known as methods of monitoring brainactivities from the outside. In the following, description will be givenmainly focusing on MEG. It is noted, however, that the present inventionis also applicable to EEG, when the magnetic field in MEG is replaced bythe electric field.

The term “electromagnetic field” originally refers to co-existence of“electric field” and “magnetic field.” In the present specification,however, an expression “observe an electromagnetic field” will be usedgenerally, where the physical amount “to be observed” is an “electricfield” or a “magnetic field.”

[First Embodiment]

FIG. 1 is a schematic illustration of an exemplary configuration of anMEG system.

An MEG 12 includes an array of multi-channel SQUID flux meter (supersensitive magnetometer), for measuring a magnetic field on the scalpsurface of a subject 10. Receiving the result of measurement by MEG 12,a computer system 20 analyses the result of measurement, and performs aprocess for estimating the position of a brain current source.

The present invention relates to software processing by computer system20. It is noted, however, that part of or all of the process may beperformed by hardware to increase the speed of processing. The softwareis not specifically limited, and it may be recorded on a recordingmedium 22 such as a CD-ROM (Compact Disk Read Only Memory) or aDVD-ROM-(Digital Versatile Disc Read Only Memory) and installed incomputer system 20, or it may be distributed through a communicationnetwork.

[Principle of Brain Current Source Estimation]

Prior to detailed description of the “brain current source estimatingmethod” in accordance with the present invention, the principle andoutline of the estimation method of the present invention will besummarized.

(Recovery of Electromagnetic Field by Current Distribution on a CurvedSurface)

It is well-known as the principle of electromagnetic shield that when acurrent source is surrounded by an ideal electromagnetic shieldingsurface formed of metal and magnetic body, electromagnetic field ceaseto exist outside the shielding surface. Specifically, by theelectromagnetic field formed by the current flowing over the shieldingsurface and a collection of small magnets constituting the magneticbody, the electromagnetic field formed by the current source outside theshielding surface is fully cancelled. Further, the small magnets can bereplaced by a virtual eddy current. When viewed from a different point,this means that an electromagnetic field identical to that formed by thecurrent source can be formed outside the shielding surface by causing anappropriate current flow over the shielding surface. On the contrary,the electric field outside the shielding surface cannot be fullyrecovered, whatever current is caused to flow over the shieldingsurface. This fact will be referred to as the “principle ofelectromagnetic field recovery.”

FIG. 2 is a schematic illustration representing a magnetic fieldgenerated by a current source, observed on an appropriate curvedsurface.

As can be seen from FIG. 2, when a virtual surface is prepared betweenan observing surface and the current source, it is possible to recoverthe electromagnetic field formed by the current source on the observingsurface by causing an appropriate flow of current on the virtualsurface, in accordance with the “principle of electromagnetic fieldrecovery.”

FIG. 3 is a schematic illustration representing a relation between acurrent source in the brain and a plurality of virtual curved surfaces.

Referring to FIG. 3, considering that the electromagnetic field formedby the current attenuates in inverse proportion to the square ofdistance, the expansion of current on the virtual curved surface(virtual curved surface 2) equivalent to the current source becomeswider as the virtual surface is further away from the current source.Therefore, the expansion of the current on virtual curved surface 1 iswider than that on virtual curved surface 2.

On a virtual curved surface (virtual curved surface 3) on the sideopposite to the observing surface with respect to the current source, itis impossible to fully recover the electromagnetic field formed on theobserving surface by the current source.

According to the present invention, based on the principle describedabove, the current source is estimated from the observed data of theelectromagnetic field on the observing surface.

(Principle of Current Source Estimation)

Assume that a magnetic field (or an electric field) formed by a currentsource generated in the brain is observed on an observing surface thatis close to the surface of the brain, as shown in FIG. 2.

A virtual curved surface in the brain is considered, and currentdistribution on the virtual curved surface that recovers the observedmagnetic field is calculated. When the virtual curved surface is movedto the inner side of the brain with the radius made gradually smaller,the expansion of the current distribution recovering the observedmagnetic becomes smaller, and it becomes the smallest when the virtualsurface encompasses the true current source. When the virtual curvedsurface is further moved to be deeper than the current source, theexpansion of current distribution comes to be wider again, and thedifference between the magnetic field generated by the current and theobserved magnetic field also becomes larger.

Accordingly, the depth of the current source can be identified byreviewing the expansion of the current distribution recovering theobserved magnetic field and the error in recovery of the magnetic field.Further, by calculating the current distribution on the virtual curvedsurface at the depth identified in this manner, the expansion of thecurrent source can also be found. The foregoing is the principle of thepresent invention.

(Identification of Current Source Depth by Bayesian Estimation)

In order to specifically implement the principle of current sourceestimation described above, in the present invention, a procedure basedon the following probability theory is employed. The procedure will besummarized in the following.

What can be observed by MEG or EEG is several to several hundreds ofmagnetic fields (electric fields) existing near the surface of thebrain. In order to approximate the current distribution on the virtualcurved surface, lattice points are set on the virtual curved surface,and a current dipole (or an appropriate current source model) isallocated to each lattice point. In order to estimate the currentdistribution with high resolution, it is necessary to increase thenumber of lattice points to increase the density of the lattice points.

When the number of lattice points on the virtual curved surface isincreased to be larger than the number of observation points, a uniquesolution cannot be determined. When estimation points is larger innumber than the observation points, the number of parameters to beestimated becomes larger than the number of observed data, and hence,the observed magnetic field would be better recovered even on a virtualsurface that is positioned deeper than the current source.

In order to cope with such problems, the current distribution on thevirtual curved surface is estimated utilizing Bayesian estimationtheory. At the time of Bayesian estimation, prior distribution thatrepresents prior information of the current source is used.Specifically, as it is considered that brain current sources existlocalized at a plurality of positions, a prior distribution that leadsto the smallest possible expansion of current distribution isintroduced. Namely, a prior distribution is introduced in which acurrent dipole on a lattice point of which magnitude cannot be wellidentified only from the observed data comes to have a magnitude closeto zero. This can be realized by introducing hierarchical priordistribution referred to as “Automatic Relevance Determination” priordistribution (Reference: R. M. Neal, Bayesian Learning for NeuralNetworks, Springer Verlag, 1996). This prior distribution will behereinafter referred to as “localized condition prior distribution.” Themanner how to select a prior distribution introducing prior informationother than the localized condition will be described later withreference to the second embodiment.

It is impossible, however, to analytically calculate posteriorprobability distribution from the localized condition prior distributionand observed data. Therefore, in the present invention, variationalBayes method is used as will be described later (Reference: H. Attias,Proc. 15th Conference on Uncertainty in Artificial Intelligence pp.21-30, 1999 and Masa-aki Sato, Neural Computation, 13, 1649-1681, 2001).It is noted that other method of approximation such as Monte Carlomethod may be used.

By Bayesian estimation using localized condition prior distribution, itbecomes possible to obtain a current distribution on the virtual curvedsurface that recovers the observed data and has smallest possibleexpansion. Further, by comparing model posterior probabilities whenestimations are made using virtual curved surfaces of different depths,the depth of the current source can be estimated.

The logarithm of the model posterior probability can be represented as asum of log-likelihood term and model complexity term. The log-likelihoodterm becomes larger as recovery error becomes smaller.

The model complexity term becomes larger when the number of effectivecurrent dipoles (that is, having a magnitude not smaller than anappropriate threshold) on the lattice points becomes smaller. As alreadydescribed, the recovery error and the expansion of the currentdistribution on the virtual curved surface (number of effective currentdipoles) become the smallest when the depth of the virtual curvedsurface matches the current source. Thus, it can be understood that themodel posterior probability becomes the largest at this time. In otherwords, the current source exists on the virtual curved surface at adepth at which the model posterior probability becomes the largest.

In summary, the present invention provides a method of estimating theposition of the brain current source, depth direction inclusive, fromthe observed data of MEG and EEG.

The basic idea of the present invention comes from the fact that anelectromagnetic field generated by a current source can be recovered bycausing an appropriate current flow over a curved surface existingbetween the current source and an observing surface, and that theexpansion of current distribution on the curved surface becomes smalleras the curved surface comes closer to the current source, as describedabove. The present invention is characterized in that, base on thisidea, the position of the current source including the depth directionis estimated by considering the fact that the model posteriorprobability becomes the largest when the curved surface encompasses thecurrent source in Bayesian estimation of the current distribution on thecurved surface that recovers the observed data, that is, by consideringthe model posterior probability.

(Where There are Plurality of Current Sources)

Though description has been made on one current source, the method isalso applicable even when there are a plurality of current sources.

The electromagnetic field generated by a current attenuates in inverseproportion to the square of distance, and therefore, the current sourceclosest to the brain surface has the largest influence on the observedmagnetic field on the brain surface. Therefore, it is possible toidentify the current sources one by one in order, starting from the oneclosest to the brain surface.

When the virtual curved surface is moved deeper from the brain surface,the model posterior probability attains the relative maximum near acurrent source closest to the brain surface (which will be referred toas the first current source). When there are two or more currentsources, there will be a plurality of local sets of current dipolescorresponding in number to the current sources, in the currentdistribution on the virtual curved surface.

The local set of current dipoles is referred to as localized currentdistribution. A local surface including individual localized currentdistribution is separated and moved in the depth direction to find themodel posterior probability. When the local surface that corresponds tothe first current source is moved, the model posterior probabilityattains to the relative maximum at the depth of the first currentsource. When other local surfaces are moved, the model posteriorprobability never attains to the relative maximum at the depth of thefirst current. Thus, the position of the first current source, that is,the position of the current source closest to the brain surface can beidentified.

In order to find the second deepest current source, the local surfacecorresponding to the first current source is fixed, the depth of theremaining local surfaces are aligned and gradually made deeper. Then themodel posterior probability attains the relative maximum at the seconddeepest current source. When the individual local surface is moved inthe depth direction again, the model posterior probability attains tothe relative maximum at the depth of the second current source, onlywhen the local surface corresponding to the second current source ismoved. In this manner, the position of the second current source can beidentified. The third and the following current sources can also beidentified in the similar manner.

The method is advantageous over the method in which the depth ofindividual local surface is moved independently, in that the time forcomputation can significantly be reduced.

(Method in Which Resolution is Increased Gradually)

According to the method of estimating brain current source of thepresent invention, it is possible to identify position of each of thecurrent sources starting from the one closest to the brain surface inthe above described manner and, in addition, it is possible to graduallyincrease the resolution of current source estimation.

First, a position of a current source is roughly estimated with a lowresolution (with a small number of lattice points on the virtual curvedsurface and a small number of sample points in the depth direction). Inthis stage, the position of a local surface corresponding to eachcurrent source is approximately determined.

Next, estimation is made with higher resolution. At this time, the areaof the local surface has become smaller as compared with the originalvirtual curved surface, and hence the resolution is naturally higherwhen the same number of lattice points are used. Accuracy in estimatingthe current source position can be improved by increasing the resolutionin the depth direction as well. If the current distribution is furtherlocalized when the resolution is made higher, it is possible to estimateagain with local surface made smaller and the resolution made stillhigher.

On the contrary, if the expansion of current distribution does not muchvary even when the resolution is improved, it means that the currentsource expands wide. In this manner, the resolution can be adjusted inaccordance with how the current source expands.

[Specific Procedure of Current Source Estimation]

In the following, specific procedure for identifying the position of a“brain current source” will be described in detail, in accordance withthe outline above.

[(I) Preparation for Current Source Estimation]

First, preparation for the process procedure for estimating the currentsource will be described.

Specifically, for the estimation of current source, the followingprocedures must be taken in advance.

(I-1) Determination of the Shape of Virtual Curved Surface and CurrentModel

(I-2) Further, in order to estimate the virtual current distributionwhile moving the virtual curved surface, it is necessary to determinesample points on the virtual curved surface and sample points in thedepth direction.

(I-3) Further, as will be described in detail later, it is necessary todetermine in advance meta parameters to designate the distribution shapeof hierarchical prior distribution for estimating current distributionas initial values.

In the following, the procedure of (I-1) Determination of the shape ofvirtual curved surface and current model will be described in graterdetail.

(I-1-1) Determination of the Shape of Virtual Curved Surface

The simplest shape of the virtual curved surface is obtained byregarding the brain as a hemisphere and assuming various hemispheres ofdifferent radii to be the virtual curved surfaces.

When a location where existence of a brain current is highly likely suchas the cerebral cortex has been known in advance by other measuringmethod such as Magnetic Resonance Imaging: MRI, the shape of the virtualcurved surface may be determined based on such information.

Particularly when the shape of the cerebral cortex is known with highprecision and it is considered that the brain current source isnon-existent in other regions, the current distribution estimation ofcerebral cortex points may be performed while omitting the search in thedepth direction, which will be described later. Further, it is alsopossible to perform the search in the depth direction only in a limitedregion.

In this case, the shapes of virtual curved surfaces having differentdepths may generally differ. It is necessary, however, to determine theshapes not to intersect with each other.

(I-1-2) Determination of Current Model

As a current model on the virtual curved surface, let us consider acurrent dipole model. It is also possible to consider other currentmodels.

Consider appropriate lattice points (sample points) {X_(n)|n=1, . . . ,N} on the virtual curved surface. Assume a current dipole of whichcurrent intensity is j_(n) on each lattice point. Here, the magneticfield formed by the current dipole j_(n) at a lattice point X_(n) on anobservation point Y_(i) (i=1, . . . , I) on the brain surface is givenby the following Biot-Savart's expression.$\mu\frac{j_{n} \times \left( {Y_{i} - X_{n}} \right)}{{{Y_{i} - X_{n}}}^{3}}$

Here, μ represents magnetic permeability, and by way of example, for avector X_(n), the expression |X_(n)| represents the absolute value ofvector X_(n). As more accurate expression, Sarvas's expression with thebrain regarded as a sphere filled with conductive solution (Reference:J. Sarvas, Phys. Med. Biol. 32, 11-22, 1987) may be used, or numericalsolution such as given by finite element method or boundary elementmethod may be used, considering detailed structure of the brain. In thefollowing, Biot-Savart's expression above will be used, for simplicityof description and representation.

Here, the magnetic field formed by the current dipole {j_(n)|n=1, . . ., N} on an observation point Y_(i) (i=1, . . . , I) is represented bythe following equation.${B\left( Y_{i} \right)} = {\sum\limits_{n = 1}^{N}{\mu\frac{j_{n} \times \left( {Y_{i} - X_{n}} \right)}{{{Y_{i} - X_{n}}}^{3}}}}$

Assuming that the direction of the magnetic field observed at theobservation point Y_(i) is a vector S_(i,c) (c=1, . . . , C), componentB_(i,c) in the direction of vector S_(i,c) of the magnetic field at thispoint can be given as$B_{i,c} = {\sum\limits_{n = 1}^{N}{\mu\frac{\left( {j_{n} \times \left( {Y_{i} - X_{n}} \right)} \right) \cdot S_{i,c}}{{{Y_{i} - X_{n}}}^{3}}}}$

Further, when a local gradient of a magnetic field is to be measured asa magnetic field to be observed, a differentiation of the equation aboveby Y_(i) will be observed.

When the direction of a current dipole at a lattice point X_(n)(position vector) is restricted, the current dipole can be representedby the following equation, with a possible independent direction of thecurrent dipole being b_(n,d)(d=1, . . . , D). In this equation, a casewhere D=3 corresponds to a case where there is no restriction on thedirection. $j_{n} = {\sum\limits_{d = 1}^{D}{j_{n,d} \cdot b_{n,d}}}$

In summary, a magnetic field formed by the current dipole at a latticepoint {X_(n)|n=1, . . . , N} of the virtual curved surface on theobservation point Y_(i) can be given as $\begin{matrix}{B_{i,c} = {\sum\limits_{n = 1}^{N}{\sum\limits_{d = 1}^{D}{j_{n,d} \cdot {G\left( {i,{c;n},d} \right)}}}}} \\{{G\left( {i,{c;n},d} \right)} = {\mu\frac{\left( {b_{n,d} \times \left( {Y_{i} - X_{n}} \right)} \right) \cdot S_{i,c}}{{{Y_{i} - X_{n}}}^{3}}}}\end{matrix}$

Here, j_(n,d) is an independent component of the current dipole at thelattice point X_(n). Further, the function G (i, c; n, d) represents thecomponent in S_(i,c) direction of the magnetic field formed by thecurrent dipole j_(n,d) at the lattice point X_(n).

The problem of estimating current distribution may be considered as aproblem of estimating a current distribution on the virtual curvedsurface {j_(n,d)|n=1, . . . , N; d=1, . . . , D} from the observedmagnetic field {B_(i,c)|i=1, . . . , I; c=1, . . . , C}.

The measured value of the magnetic field at the observation pointdescribed above may be given by the following matrix expression. Here, Gis referred to as a lead field matrix. When Sarvas's expression or anumerical solution such as obtained by the boundary element method isused in place of Biot-Savart's expression, the lead field matrix willhave different values, and other portions of the following descriptionare similarly applicable.

-   -   B=G·J        B=(B_(i,c)|i=1, . . . , I; c=1, . . . , C):(I×C) dimensional        vector        J=(J_(n,d)|n=1, . . . , N; d=1, . . . , D):(N×D) dimensional        vector        G=(G(i,c; n,d)|i=1, . . . , I; c=1, . . . , C; n=1, . . . , N;        d=1, . . . , D):(I×C)×(N×D) dimensional matrix

(Current Source Probability Model)

With the current model determined in this manner, the following “currentsource probability model” is introduced for such current distributionestimation.

(I-1-3) Current Source Probability Model

A probability model for the current model on the virtual curved surfacedescribed above will be considered.

It is assumed that the observed magnetic field is represented as a sumof the magnetic field formed by the current distribution J on thevirtual curved surface and the observation noise. Further, it is assumedthat the observation noise is Gaussian noise having an independentvariance σ² at each measuring point.

More generally, it may be possible to consider Gaussian noise having amultidimensional normal distribution, in which correlation betweennoises at respective measuring points is represented in the form of acovariance matrix. For simplicity of description, an isotropichomoscedastic noise model will be employed in the following.

Specifically, the observed magnetic field is considered as

-   -   B=G·J+ξ        J=J_(virtual curved surface)=(j_(n,d)|n=1, . . . , N; d=1, . . .        , D)        G=(G(i,c; n,d)|i=1, . . . , I; c=1, . . . , C; n=1, . . . , N;        d=1, . . . , D)        ξ=(ξ_(i,c)|i=1, . . . , I; c=1, . . . , C)        :Gaussian Noise with each Component Having Independent Variance        σ²

The probability distribution for the current model may be considered asfollows.

First, a virtual surface at a specific depth, or a set of local surfaceswith the depth of each local surface identified, will be represented byM. When a current distribution J on the virtual curved surface M isgiven, the probability P(B|J, β, M) that the observed magnetic field isB is represented by the following equation, where β=1/σ².$\begin{matrix}{{P\left( {\left. B \middle| J \right.,\beta,M} \right)} = {\exp\left\lbrack {{{- \frac{1}{2}}\beta\quad{{B - {G \cdot J}}}^{2}} + {\frac{1}{2}\left( {I \cdot C} \right)\quad\log\quad\left( {{\beta/2}\quad\pi} \right)}} \right\rbrack}} \\{\beta = {1/\sigma^{2}}}\end{matrix}$

(I-1-4) Hierarchical Prior Distribution

As already described, localized condition hierarchical priordistribution will be used as the prior distribution for the currentdistribution J on the virtual curved surface M.

The prior distribution for the current distribution J before observation(probability of J being realized) is represented by the followingequation, under the assumption of localized condition hierarchical priordistribution. $\begin{matrix}{{P_{0}\left( {\left. J \middle| \alpha \right.,\beta,M} \right)} = {\exp\left\lbrack {{{- \frac{1}{2}}\beta\quad{\sum\limits_{n = 1}^{N}{\alpha_{n}\left( {\sum\limits_{d = 1}^{D}j_{n,d}} \right)}^{2}}} + {\frac{D}{2}{\sum\limits_{n = 1}^{N}{\log\left( {\beta\quad{\alpha_{n}/2}\quad\pi} \right)}}}} \right\rbrack}} \\{{P_{0}\left( {\left. \beta \middle| \tau \right.,M} \right)} = {\Gamma\left( {\left. \beta \middle| {1/\tau} \right.,\gamma_{\beta\quad 0}} \right)}}\end{matrix}$

Here, Γ( . . . ) represents gamma distribution, which is defined below.In the expression below, Γ(γ₀) is a gamma function, which is alsodefined below. $\begin{matrix}{{\Gamma\left( {\left. \beta \middle| b \right.,\gamma_{0}} \right)} \equiv {\frac{1}{\beta}\left( {\gamma_{0}{\beta/b}} \right)^{\gamma_{0}}\frac{1}{\Gamma\left( \gamma_{0} \right)}{\mathbb{e}}^{{- \gamma_{0}}{\beta/b}}}} \\{{\Gamma\left( \gamma_{0} \right)} \equiv {\int_{0}^{\infty}{{\mathbb{d}t}\quad t^{\gamma_{0} - 1}\quad{\mathbb{e}}^{- t}}}}\end{matrix}$

Further, α and τ are hyper parameters introduced to model the currentdistribution J and prior distribution for inverse variance of noise β.Hierarchical prior distributions P₀(α|M) and P₀(τ|M) for α and τ are$\begin{matrix}{{P_{0}\left( \alpha \middle| M \right)} = {\sum\limits_{n = 1}^{N}{\Gamma\left( {\left. \alpha_{n} \middle| {\overset{\_}{\alpha}}_{0} \right.,\gamma_{\alpha\quad 0}} \right)}}} \\{{P_{0}\left( \tau \middle| M \right)} = {\Gamma\left( {\left. \tau \middle| {\overset{\_}{\tau}}_{0} \right.,\gamma_{\tau\quad 0}} \right)}}\end{matrix}$

In the equations above, meta parameters bar α₀, γ_(α0), τ₀ and barγ_(τ0) determining the distribution shape of the hierarchical priordistributions P₀(α|M) and P₀(τ|M) for α and τ will be described indetail later. Here, “bar” preceding a name of a variable means that thevariable has the sign “-” thereabove.

(I-1-5) Bayesian Estimation

When data B of a magnetic field is observed, the posterior probabilitydistribution P(J|B, M) that the current distribution is J can becalculated in the following manner, using Bayesian theorem.P(J|B,M)=∫dβ dα dτ P(J, β, α, τ|B, M)

Here, the posterior probability distribution P(J, β, α, τ|B, M) is givenas${P\left( {J,\beta,\alpha,\left. \tau \middle| B \right.,M} \right)} = \frac{{P\left( {\left. B \middle| J \right.,\beta,M} \right)}\quad{P_{0}\left( {\left. J \middle| \alpha \right.,\beta,M} \right)}\quad{P_{0}\left( {\left. \beta \middle| \tau \right.,M} \right)}\quad{P_{0}\left( \alpha \middle| M \right)}\quad{P_{0}\left( \tau \middle| M \right)}}{P\left( B \middle| M \right)}$P(B|M) = ∫𝕕J  𝕕β  𝕕α  𝕕τ  P(B|J, β, M)  P₀(J|α, β, M)P₀(β|τ, M)  P₀(α|M)  P₀(τ|M)

Using this posterior probability distribution, an expected value E[J|B,M] is given as

{overscore (J)}=E[J|B,M]=∫dJdβdαdτP(J, β, α, τ|J, M)J

Further, P(B|M) represents marginal likelihood of the virtual curvedsurface model M. When the depth of a current source is estimated, anumber of models having different depths are compared. At this time, itis assumed that there is no prior information about the depth.Specifically, in the following, it is assumed that P(M)=constant.

When the observation data B is given, the probability that the model Mis a true model, that is, the model posterior probability P(M|B), is inproportion to the model marginal likelihood P(B|M) under the assumptiondescribed above, and hence, the following relation holds.

P(M|B)∝P(B|M)

(I-1-6) Variational Bayes Method

It is generally impossible to analytically calculate the model marginallikelihood when the probability model and the hierarchical priordistribution are given.

Therefore, as a method of calculating by approximation the modelmarginal likelihood, variational Bayes method is used. It is possible touse other method of approximation, such as Monte Carlo method andLaplacian approximation.

The “variational Bayes method” will be briefly outlined, and specificprocedure will be described later.

In order to calculate a true posterior distribution P(J, β, α, τ|B, M)by approximation, a trial posterior distribution Q(J, β, α, τ) isconsidered.

Closeness between the two probability distributions P(J, β, α, τ|B, M)and Q(J, β, α, τ) can be calculated by using K-L distance given by thefollowing expressions. $\begin{matrix}{{KL}\left( {{Q\left. P \right)} = {\int\quad{{\mathbb{d}J}\quad d\quad\beta\quad d\quad\alpha\quad d\quad{{\tau Q}\left( {J,\beta,\alpha,\tau} \right)}{\log\left\lbrack \frac{Q\left( {J,\beta,\alpha,\tau} \right)}{P\left( {J,\beta,\alpha,{\tau ❘B},M} \right)} \right\rbrack}}}} \right.} \\{= {{{{\log P}\left( {B❘M} \right)} - {F(Q)}} \geq 0}} \\{{F(Q)} \equiv {\int\quad{{\mathbb{d}J}\quad d\quad\beta\quad d\quad\alpha\quad d\quad{{\tau Q}\left( {J,\beta,\alpha,\tau} \right)} \times {\log\left\lbrack \frac{{P\left( {{B❘J},\beta,M} \right)}{P_{0}\left( {{J❘\alpha},\beta,M} \right)}{P_{0}\left( {{\beta ❘\tau},M} \right)}{P_{0}\left( {\alpha ❘M} \right)}{P_{0}\left( {\tau ❘M} \right)}}{Q\left( {J,\beta,\alpha,\tau} \right)} \right\rbrack}}}}\end{matrix}$

The K-L distance attains to zero only when the two distributions areequal to each other, and otherwise it always assumes a positive value.

When a free energy F(Q) for the trial posterior distribution Q isdefined in the expression above, the following inequality results.

-   -   F(Q)≦log PB|M)

Specifically, the distribution Q(J, β, α, τ) that maximizes the freeenergy F(Q) becomes equal to the true posterior distribution P(J, β, α,τ|B, M), and at this time, the free energy is equal to the marginallog-likelihood.

According to variational Bayes method, the form of function Q(J, β, α,τ) is restricted to the following form, to maximize the free energy.

Q(J, β, α, τ)=Q_(J)(J, β)Q_(α)(α, τ)

By alternately repeating the step of fixing the second term Qα on theright side of the equation above and maximizing F(Q) with respect to thefirst term Q_(J) on the right side and the step of fixing the first termQ_(J) and maximizing F(Q) with respect to the second term Qα, adistribution Q* is obtained that attains the relative maximum of freeenergy F(Q).

[(II) Procedure of Current Source Estimation]

A specific procedure for estimating the current source after thepreparation above will be described with reference to the figures.

FIG. 4 is a flow chart representing an overall procedure of the braincurrent source estimating method in accordance with the presentinvention.

Referring to FIG. 4, first, when the process of estimating a position ofa brain current source starts (step S100), values of meta parameters fordesignating the shape of hierarchical prior distribution for estimatingthe current distribution and an initial value of a variable to beestimated are determined (step S102).

Thereafter, initial estimation of the current source is performed, usingan initial resolution, so as to extract candidates of the current source(step S104).

Then, among the current source candidates extracted in this manner, aposition of a current source that is closest to the brain surface isestimated (step S106). Thereafter, depths of other current sources isidentified successively (step S108).

After the positions of current sources are identified with the firstresolution in the above described manner, re-estimation of the positionsof the current sources is performed with the spatial resolutionincreased (step S110), and the process ends (step S112).

Processes of respective steps of FIG. 4 will be described in graterdetail in the following.

(II-1) Initial Estimation of Current Source Using Initial Resolution(Extraction of Current Source Candidates)

FIG. 5 is a flow chart representing the process of step S104 among thesteps shown in FIG. 4, that is, the process of initial estimation of thecurrent source using an initial resolution.

Referring to FIG. 5, first, sample points in the depth direction withthe initial resolution are determined to be {Rk|k=1, . . . , K}. Currentdistribution is estimated for the virtual curved surface at each depthRk. Further, based on the current distribution estimated in this manner,posterior probability for each depth Rk is calculated. The posteriorprobability corresponds to the free energy value, which will bedescribed later (step S202).

For the current distribution at the depth R_(M) at which the posteriorprobability calculated in this manner becomes the highest, a relativemaximum point of current intensity is found (step S204).

Further, a local surface that encompasses each relative maximum point isdetermined. Assuming that there are L local surfaces, each local surfacewill be the candidate of localized current source (step S206).

In the following, the process of step S102 of FIG. 4 and the process ofstep S202 of FIG. 5 that follows, will be described in grater detail.

(II-1-1) Determination of Meta Parameter Values for Designating Shape ofHierarchical Prior Distribution for Estimating Current Distribution andInitial Values of Variables to be Estimated

As described above, estimation of a current source is performed usingvariational Bayes method. Therefore, the process of step S102 shown inFIG. 4 is performed in the following manner, with the application ofvariational Bayes method.

By way of example, the meta parameters designating the shape ofhierarchical prior distribution are determined as follows.$\begin{matrix}{\gamma_{\beta 0} = 1} & {\left( {{{more}\quad{generally}},{0 \leq \gamma_{\beta 0}}} \right)} \\{\gamma_{\tau 0} = 0} & {\left( {{{more}\quad{generally}},{0 \leq \gamma_{\tau 0}}} \right)} \\{\gamma_{\alpha 0} = 0.1} & {\left( {{{more}\quad{generally}},{0 \leq \gamma_{\alpha 0}}} \right)} \\{{{\overset{\_}{\tau}}_{0} = {\kappa_{\tau} \cdot {V_{ar}(B)}}},{\kappa_{\tau} = 1}} & {\left( {{{more}\quad{generally}},{\kappa_{\tau} \sim 1}} \right)} \\{{V_{ar}(B)} = {\frac{1}{I \cdot C}{\sum\limits_{i = 1}^{I}\quad{\sum\limits_{c = 1}^{C}\quad\left( {B_{i,c} - \overset{\_}{B}} \right)^{2}}}}} & \quad \\{\overset{\_}{B} = {\frac{1}{I \cdot C}{\sum\limits_{i = 1}^{I}\quad{\sum\limits_{c = 1}^{C}\quad B_{i,c}}}}} & \quad \\{{{\overset{\_}{\alpha}}_{o} = {{\kappa_{\alpha} \cdot \frac{1}{\overset{\sim}{N}}}{{Tr}\left( {G^{\prime} \cdot G} \right)}}},{\kappa = 10}} & {\left( {{{more}\quad{generally}},{0 < \kappa_{\alpha}}} \right)} \\{\overset{\sim}{N} = {D \cdot N}} & \quad\end{matrix}$

Further, based on the equations above, each variable is initialized inthe following manner. $\begin{matrix}{\gamma_{\beta} = {{\frac{1}{2}{I \cdot C}} + \gamma_{\beta 0}}} & \quad \\{\gamma_{\tau} = {\gamma_{\tau 0} + \gamma_{\beta 0}}} & \quad \\{\gamma_{\alpha} = {\gamma_{\alpha 0} + {\frac{1}{2}D}}} & \quad \\{{\overset{\_}{\alpha}}_{n} = {\overset{\_}{\alpha}}_{0}} & \quad \\{{\overset{\_}{\tau} = {\overset{\_}{\tau}}_{0}}\quad} & \quad \\{\sum_{G}{= {G^{\prime} \cdot G}}} & \quad\end{matrix}$

(II-1-2) Specific Process of Variational Bayes Method for EstimatingCurrent Distribution

(1) Calculation of Expected Values of Parameters J, β (J-Step Process)

Here, expected values of parameters J and β are calculated. By thisprocess, the free energy F(Q) for Q_(J) is maximized.

By defining a diagonal matrix A as follows, Q_(J) is derived inaccordance with the following equations. In the following equations, anexpected value of a variable is represented by the variable name with“-” attached thereabove. $\begin{matrix}{{{\overset{\_}{A}\left( {n,{d;n^{\prime}},d^{\prime}} \right)} = \delta_{nn}},\delta_{dd},{{\overset{\_}{\alpha}}_{n}\quad\left( {n,{n^{\prime} = 1},\ldots\quad,{N;d},{d^{\prime} = 1},\ldots\quad,D} \right)}} \\{{\sum\limits^{\quad}\quad{= {\overset{\quad}{\sum_{G}}{+ \overset{\_}{A}}}}}\quad} \\{\overset{\_}{J} = {\overset{\quad}{\sum^{- 1}}\quad{\cdot G^{\prime} \cdot B}}} \\{\overset{\_}{\beta} = {\gamma_{\beta}\left\lbrack {{\frac{1}{2}{{B - {G \cdot \overset{\_}{J}}}}^{2}} + {\frac{1}{2}{\overset{\_}{J}}^{\prime}\overset{\_}{A}\quad\overset{\_}{J}} + {\gamma_{\beta 0}\overset{\_}{\tau}}} \right\rbrack}^{- 1}} \\{{Q_{J}\left( {J,\beta} \right)} = {{Q_{J}\left( {J❘\beta} \right)}{Q_{\beta}(\beta)}}} \\{{Q_{J}\left( J \middle| \beta \right)} = {\exp\left\lbrack {{{- \frac{1}{2}}{\beta\left( {J - \overset{\_}{J}} \right)}^{\prime}{\sum\left( {J - \overset{\_}{J}} \right)}} +} \right.}} \\\left. {{\frac{1}{2}\log{\sum }} + {\frac{1}{2}\overset{\sim}{N}{\log\left( {{\beta/2}\pi} \right)}}} \right\rbrack \\{{Q_{\beta}(\beta)} = {\Gamma\left( {\left. \beta \middle| \overset{\_}{\beta} \right.,\gamma_{\beta}} \right)}}\end{matrix}$

(2) Calculation of Expected Values of Hyper Parameters, that is,Calculation of Expected Values of Parameters α, τ (Process for α-Step)

Following the J step, expected values of hyper parameters α and τ arecalculated. In this step, a process is performed to maximize the freeenergy F(Q) for Q_(α).

The procedure can be represented as${\overset{\_}{\alpha}}_{n} = {{\gamma_{\alpha}\left\lbrack {{\gamma_{\alpha 0}{\overset{\_}{\alpha}}_{0}^{- 1}} + {\frac{1}{2}\overset{\_}{\beta}{\sum\limits_{d = 1}^{D}{\overset{\_}{j}}_{n,d}^{2}}} + {\sum\limits_{d = 1}^{D}{\left( \sum^{- 1} \right)\left( {n,{d;n},d} \right)}}} \right\rbrack}\quad}^{- 1}$(n = 1, …  , N)$\overset{\_}{\tau} = {\gamma_{\tau}\left\lbrack {{\gamma_{\tau 0}{\overset{\_}{\tau}}_{0}^{- 1}} + {\gamma_{\beta 0}\overset{\_}{\beta}}} \right\rbrack}^{- 1}$${Q_{\alpha}\left( {\alpha,\tau} \right)} = {{\Gamma\left( {\left. \tau \middle| \overset{\_}{\tau} \right.,\gamma_{\tau}} \right)}{\prod\limits_{n = 1}^{N}{\Gamma\left( {\left. \alpha_{n} \middle| {\overset{\_}{\alpha}}_{n} \right.,\gamma_{\alpha}} \right)}}}$

(3) Calculation of Free Energy

Using Q_(J) and Q_(α) calculated through the J step and α step asdescribed above, the free energy is calculated in the following manner.F = LP + H_(J) + H_(β) + H_(α) + H_(τ)${LP} = {{{{- \frac{1}{2}}\overset{\_}{\beta}{{B - {G \cdot \overset{\_}{J}}}}^{2}} - {\frac{1}{2}{Tr}{\sum^{- 1}{\sum_{G}{{+ \quad\frac{1}{2}}\left( {I \cdot C} \right)\quad\left( {< {\log\quad\beta} > {{- \log}\quad 2\pi}} \right)}}}}} < {\log\quad\beta}>={{\log\quad\overset{\_}{\beta}} + {\psi\left( \gamma_{\beta} \right)} - {\log\quad\gamma_{\beta}}}}$${\psi(\gamma)} \equiv {\frac{\mathbb{d}\left( {\log\quad{\Gamma(\gamma)}} \right)}{\mathbb{d}\gamma}\quad\psi\text{:}{digamma}\quad{function}}$$H_{J} = {{- {\frac{1}{2}\left\lbrack {{{Tr}{\sum^{- 1}\overset{\_}{A}}} - {\log{{\sum^{- 1}\overset{\_}{A}}}} - \overset{\sim}{N}} \right\rbrack}} - {\frac{1}{2}\overset{\_}{\beta}\quad{\overset{\_}{J}}^{\prime}\overset{\_}{A}\quad\overset{\_}{J}}}$$H_{\beta} = {{\gamma_{\beta 0}\left\lbrack {{\log\left( {\overset{\_}{\tau}\quad\overset{\_}{\beta}} \right)} - {\overset{\_}{\tau}\quad\overset{\_}{\beta}} + 1} \right\rbrack} + {\Phi\left( {\gamma_{\beta},\gamma_{\beta 0}} \right)}}$Φ(γ, γ₀) ≡ (log   Γ(γ) − γ  ψ(γ) + γ) − (log   Γ(γ₀) − γ₀log   γ₀ + γ₀)+  γ₀(ψ(γ) − log   γ)$H_{\tau} = {\gamma_{\tau 0}\left\lbrack {{\log\left( {\overset{\_}{\tau}/{\overset{\_}{\tau}}_{0}} \right)} - \left( {\overset{\_}{\tau}/{\overset{\_}{\tau}}_{0}} \right) + 1} \right\rbrack}$$H_{\alpha} = {\sum\limits_{n = 1}^{N}{\gamma_{\alpha 0}\left\lbrack {{\log\left( {{\overset{\_}{\alpha}}_{n}/{\overset{\_}{\alpha}}_{0}} \right)} - \left( {{\overset{\_}{\alpha}}_{n}/{\overset{\_}{\alpha}}_{0}} \right) + 1} \right\rbrack}}$

In this manner, the process from the J step to calculation of freeenergy described above is repeated until the value of free energy Fconverges. The value of free energy F(Q) after convergence gives anapproximation of marginal log-likelihood log P (B|M).

Further, model log-posterior probability differs only by the constantfrom log P (B|M), and therefore, the model having the maximum modelposterior probability is the same as the model of which free energy isthe highest, in accordance with the above described approximation.

By the above described procedure, posterior probability for the depth Rkcan be calculated. By performing the processes of steps S204 and S206 ofFIG. 5 accordingly, it is possible to find candidates of current sourceslocalized to respective local surfaces.

(II-2) Identification of Current Source Closest to the Brain Surface

Next, the process of step S106 of FIG. 4, that is, the process foridentifying the current source closest to the brain surface among thecandidates of current sources found in the manner as described above,will be described.

FIG. 6 is a flow chart representing a process for identifying a currentsource closest to the brain surface.

First, as the initial value, the value of a variable l is set to 1 (stepS302). Then, the variable l is compared with a possible maximum value Lof variable l (step S304), and when the variable l is not larger thanthe maximum value L, the process proceeds to step S306, and when thevariable l exceeds the maximum value L, the process proceeds to stepS324 (step S304). Specifically, the process from step S306 to step S322is repeated from l=1 to l=L.

In step S306, depth of a local surface other than the lth local surfaceis fixed at the depth R_(M) calculated in step S204 of FIG. 5.

The value of variable k is set to 1 (step S310). Thereafter, thevariable k is compared with a possible maximum value K of variable k(step S310), and when the variable k is not larger than the maximumvalue K, the process proceeds to step S312, and when the variable kexceeds the maximum value K, the process proceeds to step S320.

In step S312, the depth of the lth local surface is first determined tobe R_(k).

Thereafter, current distribution is estimated for the set of L localsurfaces (step S314).

Thereafter, the posterior probability (free energy) for this arrangementof the local surfaces is calculated (step S316). By incrementing thevalue of variable k by only 1 (step S318), the process returns to stepS310.

The process is performed for each of the local surfaces having thedepths from k=1 to k=K, and then, the depth R_(M)(l) of the first localsurface that provides the highest posterior probability is calculated(step S320). By incrementing the value of variable l only by 1 (stepS322), the process returns to step S304.

In this manner, among the depths R_(M)(l) calculated for each variablel, one closest to the brain surface (shallow) is detected, which isdenoted by l₁ (step S324).

Through the steps described above, it follows that the l₁th localsurface corresponds to the current source closest to the brain surface,the initial estimate value of the depth thereof is calculated asR_(M)(l₁), and the process terminates (step S326).

(II-3) Identification of Depth of Current Source Corresponding to eachLocal Surface

FIGS. 7 and 8 are flow charts representing a process for identifyingdepth of a current source corresponding to each local surface.

Referring to FIG. 7, the value of a variable s is set to 1 as an initialvalue (step S402). Thereafter, the variable s is compared with apossible maximum value L of variable s (step S404), and when thevariable s is not larger than the maximum value L, the process proceedsto step S406, and when the variable s exceeds the maximum value L, theprocess proceeds to step S434 (step S404). Specifically, the processfrom step S406 to step S432 is repeated from s=1 to s=L.

In step S406, when the depths of local surfaces identified so far arerepresented as {R_(m)(l₁), . . . , R_(M)(l_(s))}, the depths of theselocal surfaces are fixed at {R_(M)(l₁), . . . , R_(M)(l_(s)},respectively (step S406).

First, it is assumed that l is not equal to any of {l₁, . . . , l_(s)},and that l belongs to a set of {1, . . . , L} (step S408). Then, whetherall possible values are considered as the value of variable l or not isdetermined (step S410), and if all the possible values have beenconsidered, the process proceeds to step S430. Otherwise, the processproceeds to step S412.

In step S412, the depth of a local surface that is different from thelth local surface and not any of {l₁, . . . , l_(s))} is fixed at RM.

The value of variable k is set to 1 (step S414). Thereafter, thevariable k is compared with a possible maximum value K of variable k(step S416), and when the variable k is not larger than the maximumvalue K, the process proceeds to step S418, and when the variable kexceeds the maximum value K, the process proceeds to step S426.Specifically, the process from step S418 to step S424 is repeated fromk=1 to k=K.

In step S418, the depth of the first local surface is set to R_(k).

Thereafter, for the set of L local surfaces, current distribution isestimated (step s420).

Further, posterior probability (free energy) for this arrangement of thelocal surfaces is calculated (step 422).

Referring to FIG. 8, in step S426, after the process to k=K is finishedin the above described manner, the depth R_(M)(l) of the lth localsurface attaining the highest posterior probability is calculated (stepS426).

Then, other variable l that is not equal to any of {l₁, . . . , l_(s)}belongs to the set of {1, . . . , L} and is not yet processed isselected (step S428), and the process returns to step S410.

Through the above described steps, among R_(M)(l) values correspondingto all the processed variables l, the value 1 that is closest to thebrain surface is calculated and denoted by l_(s+1) (step S430).

Further, s is replaced by s+1 (step S432), and the process returns tostep S404.

When the process ends for all the possible values s, identification ofthe depths of current sources corresponding to respective local surfacesis finished (step S434).

(II-4) Re-Estimation with Increased Resolution

FIGS. 9 and 10 are flow charts representing the process of step S110shown in FIG. 4, that is, the process for re-estimating a position of acurrent source with higher spatial resolution.

Referring to FIGS. 9 and 10, first, numbers of local surfacescorresponding to the current sources estimated using the initialresolution are re-numbered so that the surfaces have the numbers 1, 2,3, . . . , L starting from the one closest to the brain surface (stepS502). Then, corresponding to the expansion of current distribution oneach local surface, the local surface is made smaller (step S504).

The depth of the local surface R_(M)(l) calculated by the process up tostep S108 of FIG. 4 is denoted as R_(M) ^(old)(l) (step S506).

New resolution and search width in the depth direction are respectivelyrepresented as ΔR and (K_(L)·ΔR). Further, resolution of lattice pointson the local surface is also increased (step S508).

The value of a variable l is set to 1 as an initial value (step S510).Thereafter, the variable l is compared with a possible maximum value Lof variable l (step S512), and when the variable l is not larger thanthe maximum value L, the process proceeds to step S514, and when thevariable l exceeds the maximum value L, the process proceeds to stepS534 (step S512). Specifically, the process from step S514 to step S532is repeated from l=1 to l=L.

In step S514, the depth of a local surface l′ other than the lth localsurface is fixed to R_(M) ^(old)(l′).

Further, the process from step S518 to step S526 is performed with k=0,±1, . . . , ±K_(L).

First, the value of variable k is set to 0 (step S516). Thereafter, theabsolute value of variable k is compared with the possible maximumabsolute value of K_(L) (step S518), and when the absolute value ofvariable k is not larger than the maximum value K_(L), the processproceeds to step S520, and when the absolute value of variable k exceedsthe maximum value K_(L), the process proceeds to step S528 (step S518).

In step S520, the depth of the lth local surface is set to (R_(M)^(old)(l)+k·ΔR).

Thereafter, for the set of L local surfaces, current distribution isestimated using the present resolution (step S522).

Further, posterior probability (free energy) for this arrangement of thelocal surfaces is calculated (step S524). Then, the value k is set tothe next one of {±1, . . . , ±K_(L)}, and the process returns to stepS518.

After the above described process is performed until k=±K_(L), the valuek that provides the highest posterior probability is calculated in stepS528, which value is denoted by k_(M).

Then, R_(M) ^(old)(l) is replaced by R_(M) ^(old)(l)+k_(M)·ΔR (stepS530). The value of variable l is incremented by only 1 (step S532), andthe process returns to step S512.

When the above described process has been performed until the value ofvariable l attains to the maximum value L and the resolution has reachedthe final resolution (step S534), the process is terminated (step S538).

When the resolution is not yet the final resolution (step S534), theresolution of the lattice points on the local surface and the resolutionin the depth direction are increased (step S536), and the processreturns to step S510.

By the “method of estimating brain current source” as described above,it becomes possible to estimate the position of a brain current source,including the position in the depth direction, from observation data ofMEG (or EEG). Further, such estimation in the depth direction isapplicable even when there are a plurality of current sources. Stillfurther, the method is applicable when the current sources are localizedas in the case of current dipoles or when the current source has anexpansion. In addition, by the method, it is possible to estimate howthe current source expands.

By additionally performing the process described with reference to FIGS.9 and 10 after the estimation with the initial resolution, it becomespossible to successively increase resolution for estimating a position.This also means that it is possible to review with the scope of searchrestricted based on physiological findings or data obtained by otherobservation method.

[Result of Simulation]

In the following, simulation results will be described, in which currentsource positions of an assumed model were estimated in accordance withthe method of estimating brain current source described above.

FIG. 11 is a top view of a magnetic field distribution in the radialdirection observed on a surface of a hemisphere, assuming that the humanbrain is a hemisphere having the radius of 10.0 (arbitrary unit).

FIG. 11 shows a magnetic field distribution on a surface of thehemisphere in which a single current source is positioned at a radius of7.0 from the center, as will be described later. In the simulation thatwill be discussed below, noise having the S/N ratio of 0.1 is added tothe observation data of the magnetic field.

FIG. 12 represents results of initial estimation of current sourcesusing initial resolution.

FIG. 12 shows the result of calculation of current distribution, inwhich, as described with reference to the process of steps S102 to S104of FIG. 4, in order to perform initial estimation of the current sourceusing the initial resolution, the depth (radius) is changed stepwise andon the hemisphere of each depth, the current distribution that attainsthe maximum free energy is calculated in accordance with variationalBayes method.

FIG. 12(a) shows the result where the radius is R=5.0, 12(b) shows theresult where the radius is R=6.0 and 12(c) shows the result where theradius is R=7.0.

FIG. 13 represents results of initial estimation of current sourcesusing initial resolution, where the radius is larger (that is, closer tothe surface of the brain).

FIG. 13(d) shows the result where the radius is R=7.5, 13(e) shows theresult where the radius is R=8.0 and 13(f) shows the result where theradius is R=9.0.

FIG. 14 represents magnitude of free energy on each virtual hemisphere,obtained when current distribution that attains maximum free energy iscalculated for virtual hemispheres at various depths assumed in thebrain.

Referring to FIG. 14, it can be understood that the maximum point offree energy exists between the radii of 7.0 and 8.0, when calculation ismade assuming that the hemisphere as a whole is a virtual curvedsurface.

As can be seen from FIG. 13(d), one point of relative maximum exists inthe current distribution on the virtual hemispherical surface having theradius of R=7.5.

As the current source is initially estimated using the initialresolution in this manner, on the virtual curved surface having themaximum free energy, a local surface that encompasses the relativemaximum point of current distribution is calculated, and the process foridentifying the current source is performed on the local surface.

FIGS. 15 and 16 represent processes for identifying the current sourceexecuted further on a local surface including a point of relativemaximum, with the resolution of lattice points on the local surface andthe resolution in the depth direction increased.

FIGS. 15(a), 15(b) and 15(c) represent current distribution on localsurfaces where the radius is R=5.0, R=6.0 and R=7.0, respectively.Similarly, FIGS. 16(d), 16(e) and 16(f) represent current distributionon local surfaces where the radius is R=7.5, R=8.0 and R=9.0,respectively.

FIG. 17 represents the free energy calculated with the depth of localsurface varied.

Referring to FIG. 17, the free energy as a result of calculation usingthe local surface has the relative maximum value near the radius of R=7,and from the result, it is possible to identify that one current sourceexists at the position of R=7.

The current distribution of the current source at this time is as shownin FIG. 15(c), and as already described, it can be understood that thecurrent distribution also has the narrowest expansion.

FIG. 18 is a top view of a magnetic field distribution observed on asurface of the brain assumed as a hemisphere, when there are two currentsources in the brain.

In the example shown in FIG. 18, it is assumed that the current sourcesare positioned at the radii of R=6.0 and R=8.0.

FIGS. 19 and 20 represent results of initial estimation of currentsources using initial resolution.

Specifically, FIGS. 19 and 20 illustrates the process (step S104) forperforming initial estimation of current sources using the initialresolution shown in FIG. 4, when there are two current sources in thebrain.

FIGS. 19(a), 19(b) and 19(c) represent current distribution on virtualhemispherical surfaces where the radius is R=5.0, R=6.0 and R=7.0,respectively.

FIGS. 20(d), 20(e) and 20(f) represent current distribution on virtualhemispherical surfaces where the radius is R=7.5, R=8.0 and R=9.0,respectively.

FIG. 21 represents radius dependency of free energy, when currentdistribution that attains maximum free energy is calculated for eachvirtual hemispherical surface of the whole hemispherical surfaces.

As shown in FIG. 21, when the free energy is calculated over the entirehemispherical surface, it is understood that the maximum value of freeenergy exists between the radii of R=7 and R=8, that is, near the radiusof R=7.5.

As can be seen from FIG. 20(d), on the virtual hemispherical surfacehaving the radius of R=7.5, there are two relative maximum points in thecurrent distribution and two local surfaces are determined.

Next, the result of calculation corresponding to the process of stepS1106 shown in FIG. 4 will be described.

FIGS. 22 and 23 illustrates the process for calculating the depth of alth local surface R_(M)(l) attaining maximum posterior probability, toidentify a current source closest to the brain surface.

Therefore, FIGS. 22 and 23 represent current distribution when a localsurface corresponding to one current source is moved while the otherlocal surface is fixed on a surface that attains the maximum free energycalculated over the entire spherical surface.

FIGS. 22(a), 22(b), 22(c) represent current distribution on localsurfaces where the radius is R=5.0, R=6.0 and R=7.0, respectively.

FIGS. 23(d), 23(e) and 23(f) represent current distribution on localsurfaces where the radius is R=7.5, R=8.0 and R=9.0, respectively.

FIG. 24 represents local surface position (radius) dependency of freeenergy when a virtual local surface is moved as shown in FIGS. 22 and23.

As shown in FIG. 24, when a local surface corresponding to one currentsource is moved to the position of R=8 while the other local surface isfixed on a surface that attains the maximum free energy calculated overthe entire spherical surface, the free energy is maximized.

Therefore, from the result shown in FIG. 24, it can be seen that theposition of the current source closest to the brain surface is at theradius of R=8.0.

Thereafter, the depth of the other current source is identified.

Here, the depth of the local surface corresponding to one current sourceis fixed at the radius of R=8.0, and the depth of the local surfacecorresponding to the other current source is moved.

FIGS. 25 and 26 represent current distribution when one local surface isfixed on a specific surface and a local surface corresponding to theother current source is moved.

FIGS. 25(a), 25(b) and 25(c) represent current distribution on localsurfaces where the radius is R=5.0, R=5.5 and R=6.0.

FIGS. 26(d), 26(e) and 26(f) represent current distribution on localsurfaces where the radius is R=6.5, R=7.0 and R=7.5.

FIG. 27 represents local surface position (radius) dependency of freeenergy when a virtual local surface is moved as shown in FIGS. 25 and26.

As can be seen from FIG. 27, the free energy is maximized at theposition of R=6.

It can be understood that by such a process, even when there are twocurrent sources in the brain, it is possible to identify the depth ofeach of the current sources.

[Second Embodiment]

In the following, the second embodiment of the present invention will bedescribed, in which, among the processes implemented by software (orpartially by hardware) of the “brain current source estimating method,”“brain current source estimating program,” and “brain current sourceestimating apparatus” in accordance with the first embodiment of thepresent invention described above, the processes from “(I) Preparationfor current source estimation” to “(3) Calculation of free energy” aremodified.

As a background for describing the modifications in accordance with thesecond embodiment, the concept of the first embodiment will besummarized in the following, and of the basic concept, portions that aremodified in the second embodiment will be described.

(Concept of Virtual Curved Surface Estimation in the First Embodiment)

In the first embodiment, as a current model on the virtual curvedsurface, current dipoles on lattice points are assumed, and a process ofsolving an inverse problem that the current distribution J is estimatedfrom observed magnetic field B on the observation surface is performed.Here, such an inverse problem is a so called “ill-posed problem” and itis difficult to find an exact solution, as (the number of lattice pointsat which current dipoles exist) is generally larger than (the number ofobservation points).

In the first embodiment, in order to solve such a problem, the method of“Bayesian estimation” has been used.

Bayesian estimation utilizes the fact that posterior probabilitydistribution when virtual surface M and observed magnetic field B aregiven can be represented by the following equation.${P\left( {\left. J \middle| B \right.,M} \right)} = \frac{{P\left( {\left. B \middle| J \right.,M} \right)}{P\left( J \middle| M \right)}}{P\left( B \middle| M \right)}$

Here, P(J|B,M) is the probability distribution on the virtual curvedsurface M under the condition that the magnetic field B has beenobserved, and it represents the “posterior probability distribution.”P(B|J, M) is “data likelihood,” P(J|M) is prior distribution and P(B|M)is marginal likelihood.

The marginal likelihood is given by

P(B|M)=∫dJP(B|J,M)P(J|M)

The posterior probability P(M|B) is given by${P\left( M \middle| B \right)} = {\frac{{P\left( B \middle| M \right)}{P(M)}}{P(B)} \propto {P\left( B \middle| M \right)}}$

Assuming that there is no prior information related to the depth,P(M)=constant, and therefore, the probability that the model M is thetrue model when the observation data B is given is in proportion to themodel marginal likelihood P(M|B) under the assumption above, as alreadydescribed with reference to the first embodiment.

Specifically, in the first embodiment, when the current model isselected, a model that attains the highest posterior probability, or thehighest marginal likelihood, has been selected.

For such a process, the marginal log-likelihood log P(B|M) representedby the following equation is the object of calculation.log P(B|M)=<log(P(B|J,M))>_(J) −KL(P(J|B,M)∥P(J|M))=(expectedlog-likelihood)−(distance of prior/posterior distribution)˜−(recoveryerror)−(effective degree of freedom of parameters)

As represented by the equation above, that the marginal log-likelihoodis maximized corresponds to a concept that the error and the effectivedegree of freedom of parameters are minimized, in other words, that theerror and the expansion of current distribution are minimized.

In the first embodiment, the problem of maximizing the marginallog-likelihood described above has been solved by the so called“variational Bayes method” in which the trial posterior distribution Qis introduced as an approximation of the posterior distribution, and theposterior distribution is determined to maximize the free energy (O).

(Concept of Prior Distribution in the First Embodiment)

In the first embodiment, as a prior information of the brain currentsource, a localized condition that the current sources exist at aplurality of positions in the brain is used, which is represented in theform of hierarchical prior distribution. Specifically, hyper parametersα_(n)(n=1, . . . , N) are introduced and the form of the hierarchicalprior distribution given by the following expressions are assumed.${P_{0}\left( {\left. J \middle| \alpha \right.,\beta} \right)} \propto {\exp\left\lbrack {{- \frac{1}{2}}\beta{\sum\limits_{n = 1}^{N}{\alpha_{n}{J_{n}}^{2}}}} \right\rbrack}$${P_{0}\left( {\left. \alpha_{n} \middle| {\overset{\_}{\alpha}}_{0} \right.,\gamma_{\alpha 0}} \right)} \propto {\exp\left\lbrack {{{- \gamma_{\alpha 0}}{\alpha_{n}/{\overset{\_}{\alpha}}_{0}}} + {\left( {\gamma_{\alpha 0} - 1} \right)\log\quad\alpha_{n}}} \right\rbrack}$

From the expressions above, it can be understood that the hyperparameter an is in proportion to inverse variance (the reciprocal ofvariance) of current J_(n). When the values of variance of current J_(n)and inverse variance β of noise are known in advance, a priordistribution may be used in which the value α_(n) is fixed to the valuecalculated from the values of variance of current and inverse varianceof noise. The value of variance of the current, however, is generallyunknown, and therefore, in the first embodiment, a hierarchical priordistribution P₀(α_(n)|bar α₀ γ_(α0)) for the hyper parameter α_(n) isassumed, and Bayesian estimation is performed on α_(n).

Here, bar α_(n) and γ_(α0) are meta parameters for determining the shapeof hierarchical prior distribution, and common values are assumed forevery α_(n). When Bayesian estimation is performed using the localizedcondition hierarchical prior distribution, Bayesian estimation thatrecovers the observed data with smallest possible number of currentdipoles is performed as described above.

Minimum norm estimation method, which is often used in current sourceestimation of MEG (Reference: M. S. Hämäläinen and R. J. Ilmoniem, Med.& Biol. Eng. & Comput., 32, 35-42, 1994) corresponds, from the viewpoint of Bayesian estimation, to a process in which the value α_(n) asthe prior distribution is made common to all the points, and the valueis directly determined to be an appropriate value.

In the minimum norm estimation, total sum of current intensities at allthe points is to be minimized, and hence points having small currentintensities appear in large number in the estimated solution. Therefore,this method is disadvantageous in that it is difficult to determinewhether such small current distributions represent the true currentdistribution or noise derived from estimation error.

The localized condition hierarchical prior distribution solves thisproblem of the minimum norm estimation method.

(Modifications in the Second Embodiment)

In the second embodiment, prior distribution considering not only thelocal condition but also continuous condition is assumed. When currentsource distribution is to be estimated with high precision with thespatial resolution in the order of millimeters, it is necessary toconsider continuity of current distribution. Actually, neural cells inthe cerebral cortex has a columnar structure having the radius of about10 mm, and when a neural activity takes place, it is likely that neuralcells in an area having the radius of about 20 mm fire simultaneously.Considering such points, in the second embodiment, a hierarchical priordistribution is assumed that combines the local condition and continuouscondition, as will be described in the following.

In order to introduce the continuous condition, an internal variable Zis newly introduced, and the current J is represented as the internalvariable Z smoothed through a smoothing filter W. Here, though notlimiting, a Gauss filter may be used as such a smoothing filter W. Now,the prior distribution can be represented as${P_{0}\left( {\left. J \middle| Z \right.,\alpha,\beta} \right)} \propto {\exp\left\lbrack {{- \frac{1}{2}}\beta{\sum\limits_{n = 1}^{N}{\alpha_{n}\left( {J_{n} - {W \cdot Z_{n}}} \right)}^{2}}} \right\rbrack}$${P_{0}\left( Z \middle| \lambda \right)} \propto {\exp\left\lbrack {{- \frac{1}{2}}\beta{\sum\limits_{n = 1}^{N}{\lambda_{n}{Z_{n}}^{2}}}} \right\rbrack}$

In the prior distribution above, integration of Z can be executedeasily, as represented by${P_{0}\left( {{J❘\alpha},\beta,\lambda} \right)} = {{\int{{\mathbb{d}Z}\quad{P_{0}\left( {{J❘Z},\alpha,\beta} \right)}{P_{0}\left( {Z❘\lambda} \right)}}} \propto {\exp\left\lbrack {{- \frac{1}{2}}\beta\quad{J^{\prime}\left( {A^{- 1} + {W \cdot \Lambda^{- 1} \cdot W^{\prime}}} \right)}^{- 1}J} \right\rbrack}}$where A and Λ represent diagonal matrixes having {α_(n)|n=1, . . . , N}and {λ_(n)|n=1, . . . , N} as diagonal components, respectively.

Specifically, assuming the prior distribution P₀(J|Z, α, β)P₀(Z|λ) byintroducing the internal variable Z is the same as assuming the priordistribution that has a covariance matrix β⁻¹(A⁻¹+W·Λ⁻¹·W′) for thecurrent J. It is easier, however, to apply variational Bayes method whenthe prior distribution is represented by using the internal variable Z.Therefore, the prior distribution will be represented by using theinternal variable Z in the following. Further, the value of hyperparameter λ_(n) is not known in advance as in the case of hyperparameter α_(n), Bayesian estimation is performed onα_(n) and λ_(n),assuming the hierarchical prior distribution of the following form.

P₀(α_(n)|{overscore (α)}_(0n), γ_(α0n))∝exp [−γ_(α0n)α_(n)/{overscore(α)}_(0n)+(γ_(α0n)−1)log α_(n)]

P₀(λ_(n)|{overscore (λ)}_(0n), γ_(λ0n))∝exp [−γ_(λ0n)λ_(n)/{overscore(λ)}_(0n)+(γ_(λ0n)−1)log λ_(n)]

The hierarchical prior distribution extended in this manner encompassesthe localized condition prior distribution of the first embodiment.Specifically, when the smoothing filter matrix W is identically 0, thehierarchical prior distribution combining the continuous condition andthe localized condition coincides with the localized condition priordistribution of the first embodiment.

It is assumed that meta parameters bar α_(n0), γ_(α0n), bar λ_(0n) andγ_(λ0n) for determining the shape of hierarchical prior distribution maydepend more generally on the positions of lattice points. When there isonly the observation data of MEG (or EEG), there is no prior informationrelated to these values. Therefore, as in the first embodiment, theseparameters are substantially determined commonly without any dependenceon the location of lattice points, as will be described in thefollowing. When observation data obtained by other observation means isavailable, however, the values of meta parameters may be determined foreach lattice point using such information, as will be described withreference to the third embodiment.

(Specific Calculation Procedure of the Second Embodiment)

The calculation procedure described in the following is in most partsimilar to that of the first embodiment. Therefore, in the following,description will be made mainly focusing on the modifications from thefirst embodiment.

First, as in the first embodiment, when a current distribution J on avirtual curved surface is given and the observed magnetic field is B, amodel probability distribution P(B|J, β) is given, using logarithmicrepresentation, by the following equation.${\log\quad{P\left( {{B❘J},\beta} \right)}} = {{{- \frac{1}{2}}\beta{{B - {G \cdot J}}}^{2}} + {\frac{1}{2}\left( {I \cdot C} \right){\log\left( {{\beta/2}\pi} \right)}}}$

-   observed magnetic field B: (I×C) dimensional vector-   current distribution J: R dimensional vector, R=N×D, N: number of    lattice points, D: number of components-   lead field matrix: G(I×C)×R dimensional matrix

[Hierarchical Prior Distribution]

The hierarchical prior distribution in accordance with the secondembodiment is represented by the following equations.${\log\quad{P_{0}\left( {{J❘Z},\beta,\alpha} \right)}} = {{{- \frac{1}{2}}{\beta\left( {J - {W\quad Z}} \right)}^{\prime}{A\left( {J - {W\quad Z}} \right)}} + {\frac{D}{2}{\sum\limits_{n = 1}^{N}{\log\left( {\beta\quad{\alpha_{n}/2}\pi} \right)}}}}$${\log\quad{P\left( {{Z❘\beta},\lambda} \right)}} = {{{- \frac{1}{2}}\beta\quad Z^{\prime}\Lambda\quad Z} + {\frac{D}{2}{\sum\limits_{n = 1}^{N}{\log\left( {\beta\quad{\lambda_{n}/2}\pi} \right)}}}}$(A)(n, d; n^(′), d^(′)) = δ_(nn^(′))δ_(dd^(′))α_(n)(Λ)(n, d; n^(′), d^(′)) = δ_(nn^(′))δ_(dd^(′))λ_(n)(n, n^(′) : 1, …  , N; d, d^(′) : 1, …  , D)${\log\quad{P_{0}(\alpha)}} = {\sum\limits_{n = 1}^{N}\left\lbrack {{{- \gamma_{\alpha\quad 0\quad n}}{\overset{\_}{\alpha}}_{0\quad n}^{- 1}\alpha_{n}} + {\left( {\gamma_{\alpha\quad 0\quad n} - 1} \right)\log\quad\alpha_{n}} - {\Phi_{G}\left( {\gamma_{\alpha\quad 0\quad n},{\overset{\_}{\alpha}}_{0\quad n}^{- 1}} \right)}} \right\rbrack}$${\log\quad{P_{0}(\lambda)}} = {\sum\limits_{n = 1}^{N}\left\lbrack {{{{- \gamma_{\lambda\quad 0\quad n}}{\overset{\_}{\lambda}}_{0\quad n}^{- 1}\lambda_{n}} + {\left( {\gamma_{{\lambda 0}\quad n} - 1} \right)\log\quad\lambda_{n}} - {\Phi_{G}\text{(}\gamma_{\lambda\quad 0\quad n}}},{{\overset{\_}{\lambda}}_{0\quad n}^{- 1}\text{)}}} \right\rbrack}$log   P₀(β❘τ) = −γ_(β  0)τβ + (γ_(β0) − 1)log   β − Φ_(G)(γ_(β0), τ)${\log\quad{P_{0}(\tau)}} = {{{- \gamma_{\tau 0}}{\overset{\_}{\tau}}_{0}^{- 1}\tau} + {\left( {\gamma_{\tau 0} - 1} \right)\log\quad\tau} - {\Phi_{G}\left( {\gamma_{\tau 0},{\overset{\_}{\tau}}_{0}^{- 1}} \right)}}$Φ_(G)(γ, τ) = log   Γ(γ) − γlog(γ  τ)

Here, as in the first embodiment, β represents inverse variance ofnoise.

[Calculation of Trial Posterior Distribution]

At this time, it is assumed that the trial posterior distribution isrepresented by the form of a product of Q_(J), Q_(Z) and Q_(α), such asQ=Q_(J)(J, β)Q_(Z)(Z)Q_(α)(α, λ, τ).

Then, the free energy F(Q) is maximized successively for each of Q_(J),Q_(Z) and Q_(α). First, in the first step, Q_(Z) and Q_(α) are fixed,and F(Q) is maximized for Q_(J), and the following equations result.Q_(J)(J, β) = Q_(J)(J❘β)Q_(β)(β)${\log\quad{Q_{J}\left( {J❘\beta} \right)}} = {{{- \frac{1}{2}}{\beta\left( {J - \overset{\_}{J}} \right)}^{\prime}{\Sigma\left( {J - \overset{\_}{J}} \right)}} + {\frac{1}{2}\log{\Sigma }} + {\frac{1}{2}\left( {I \cdot C} \right){\log\left( {{\beta/2}\pi} \right)}}}$${{\log\quad{Q_{\beta}(\beta)}} = {{{- \gamma_{\beta}}{\overset{\_}{\beta}}^{- 1}\beta} + {\left( {\gamma_{\beta} - 1} \right)\log\quad\beta} - {\Phi_{G}\text{(}\gamma_{\beta}}}},{{\overset{\_}{\beta}}^{- 1}\text{)}}$

Next, in the second step, Q_(J) and Q_(α) are fixed, and F(Q) ismaximized for Q_(Z), and the following equation results.${\log\quad{Q_{Z}(Z)}} = {{{- \frac{1}{2}}{\overset{\_}{\beta}\left( {Z - \overset{\_}{Z}} \right)}^{\prime}{\sum\limits_{Z}\left( {Z - \overset{\_}{Z}} \right)}} + {\frac{1}{2}\log{\sum\limits_{Z}}} + {\frac{1}{2}R\quad{\log\left( {{\overset{\_}{\beta}/2}\pi} \right)}}}$

Further, in the third step, Q_(J) and Q_(Z) are fixed, and F(Q) ismaximized for Q_(α), and the following equations result.Q_(α)(α, λ, τ) = Q_(α)(α)Q_(λ)(λ)Q_(τ)(τ)${\log\quad{Q_{\alpha}(\alpha)}} = {\sum\limits_{n = 1}^{N}\left\lbrack {{{- \gamma_{\alpha\quad n}}{\overset{\_}{\alpha}}_{n}^{- 1}\alpha_{n}} + {\left( {\gamma_{\alpha\quad n} - 1} \right)\log\quad\alpha_{n}} - {\Phi_{G}\left( {\gamma_{\alpha\quad n},{\overset{\_}{\alpha}}_{n}^{- 1}} \right)}} \right\rbrack}$${\log\quad{Q_{\lambda}(\lambda)}} = {\sum\limits_{n = 1}^{N}\left\lbrack {{{{- \gamma_{\lambda\quad n}}{\overset{\_}{\lambda}}_{n}^{- 1}\lambda_{n}} + {\left( {\gamma_{\lambda\quad n} - 1} \right)\log\quad\lambda_{n}} - {\Phi_{G}\text{(}\gamma_{\lambda\quad n}}},{{\overset{\_}{\lambda}}_{n}^{- 1}\text{)}}} \right\rbrack}$${\log\quad{Q_{\tau}(\tau)}} = {{{- \gamma_{\tau}}{\overset{\_}{\tau}}^{- 1}\tau} + {\left( {\gamma_{\tau} - 1} \right)\log\quad\tau} - {\Phi_{G}\left( {\gamma_{\tau},{\overset{\_}{\tau}}^{- 1}} \right)}}$

The specific method of calculating unknown numbers in the equationsabove will be described later.

As in the first embodiment, a trial posterior distribution that attainsthe relative maximum of free energy F(Q) is found again throughrepetitive calculation in the second embodiment.

Specifically, by the procedure described below, the trial posteriordistribution Q is calculated through the first to third steps, and thisprocedure is repeated until the value of free energy F(Q) converges. Howto determine the constant term at this time will be described later.

(1) Procedure of the First Step (Referred to as J-Step)

Expected values of current J and inverse variance β, and variable γ_(β)are calculated in accordance with the following equations.$\Sigma = {{{G^{\prime}G} + \overset{\_}{A}} = {\sum\limits_{G}{+ \overset{\_}{A}}}}$$\overset{\_}{J} = {\sum\limits^{- 1}\left( {{G^{\prime} \cdot B} + {\overset{\_}{A}\quad W\quad\overset{\_}{Z}}} \right)}$γ_(β) = (I ⋅ C)/2 + R/2 + γ_(β  0)${\gamma_{\beta}{\overset{\_}{\beta}}^{- 1}} = {{\frac{1}{2}{{B - {G \cdot \overset{\_}{J}}}}^{2}} + {\frac{1}{2}\left( {\overset{\_}{J} - {W\quad\overset{\_}{Z}}} \right)^{\prime}{\overset{\_}{A}\left( {\overset{\_}{J} - \quad{W\quad\overset{\_}{Z}}} \right)}} + {\frac{1}{2}{\overset{\_}{Z}}^{\prime}\overset{\_}{\Lambda}\quad\overset{\_}{Z}} + {\frac{R}{2}{\overset{\_}{\beta}}^{- 1}} + {\gamma_{\beta 0}\overset{\_}{\tau}}}$

(2) Procedure of the Second Step (Z-Step)

Thereafter, expected value of internal variable Z is calculated asfollows.Σ_(Z)={overscore (Λ)}+W′{overscore (A)}W{circumflex over (Σ)}_(Z) =I+({overscore (Λ)}⁻¹ W′{overscore (Λ)}) W{overscore (Z)}=Σ _(Z) ⁻¹ W′{overscore (A)}{overscore (J)}

(3) Procedure of the Third Step (α-Step)

Thereafter, the following calculations are made for the hyper parametersα_(n), γ_(αn), λ_(n), γ_(λn), γ_(τ).$\gamma_{\alpha\quad n} = {\frac{D}{2} + \gamma_{{\alpha 0}\quad n}}$${\gamma_{\alpha\quad n}{\overset{\_}{\alpha}}_{n}^{- 1}} = {{\frac{1}{2}\overset{\_}{\beta}{{{\overset{\_}{J}}_{n} - \left( {W\quad\overset{\_}{Z}} \right)_{n}}}^{2}} + {\frac{1}{2}{\sum\limits_{d = 1}^{D}\left\lbrack {{\left( \sum\limits^{- 1} \right)\left( {{n.d};{n.d}} \right)} + {\left( {W{\sum\limits_{Z}^{- 1}W^{\prime}}} \right)\left( {n,{d;n},d} \right)}} \right\rbrack}} + {\gamma_{{\alpha 0}\quad n}{\overset{\_}{\alpha}}_{0\quad n}^{- 1}}}$$\gamma_{\lambda\quad n} = {\frac{D}{2} + \gamma_{{\lambda 0}\quad n}}$${\gamma_{\lambda\quad n}{\overset{\_}{\lambda}}_{n}^{- 1}} = {{\frac{1}{2}\overset{\_}{\beta}{{\overset{\_}{Z}}_{n}}^{2}} + {\frac{1}{2}{\sum\limits_{d = 1}^{D}{\left( \sum\limits_{Z}^{- 1} \right)\left( {{n.d};{n.d}} \right)}}} + {\gamma_{{\lambda 0}\quad n}{\overset{\_}{\lambda}}_{0\quad n}^{- 1}}}$γ_(τ) = γ_(β0) + γ_(τ  0)${\gamma_{\tau}{\overset{\_}{\tau}}^{- 1}} = {{\gamma_{\beta 0}\overset{\_}{\beta}} + {\gamma_{\tau 0}{\overset{\_}{\tau}}_{0}^{- 1}}}$

The trial posterior distribution Q can be calculated through theprocedures (1) to (3) above.

[Calculation of Free Energy]

The free energy of the second embodiment is calculated in accordancewith the following equation.F=LP+H _(J) +H _(Z) +H _(α) +H _(λ) +H _(β) +H _(τ)

Expected log-likelihood LP and model complexity terms H_(J), H_(Z),H_(α), H_(λ), H_(β), H_(τ) are calculated as follows.${LP} = {{{- \frac{1}{2}}\overset{\_}{\beta}{{B - {G \cdot \overset{\_}{J}}}}^{2}} - {\frac{1}{2}{{Tr}\left( {\sum\limits^{- 1}{G^{\prime}G}} \right)}} + {\frac{1}{2}\left( {I \cdot C} \right)\left( {{\text{<}\log\quad\beta\text{>}} - {\log\quad 2\pi}} \right)}}$${\text{<}\log\quad\beta\text{>} = \log\quad\overset{\_}{\beta}} + {\psi\left( \gamma_{\beta} \right)} - {\log\quad\gamma_{\beta}}$$H_{J} = {{\frac{1}{2}\left\lbrack {{\log{{\overset{\_}{A}\sum\limits^{- 1}}}} - {{Tr}\left( {\overset{\_}{A}\sum\limits^{- 1}} \right)} + R} \right\rbrack} - {\frac{1}{2}{\overset{\_}{\beta}\left( {\overset{\_}{J} - \quad{W\quad\overset{\_}{Z}}} \right)}^{\prime}{\overset{\_}{A}\left( {\overset{\_}{J} - {W\quad\overset{\_}{Z}}} \right)}} - {\frac{1}{2}{Tr}{\sum\limits_{Z}^{- 1}{W^{\prime}\overset{\_}{A\quad}W}}}}$$H_{\beta} = {{\gamma_{\beta 0}\left\lbrack {{\log\left( {\overset{\_}{\beta}\overset{\_}{\tau}} \right)} - {\overset{\_}{\beta}\overset{\_}{\tau}} + 1} \right\rbrack} + {\Phi\left( {\gamma_{\beta},\gamma_{\beta 0}} \right)}}$$H_{Z} = {{\frac{1}{2}\left\lbrack {{\log{\overset{- 1}{\hat{\sum\limits_{Z}}}}} - {{Tr}\left( \overset{- 1}{\hat{\sum\limits_{Z}}} \right)} + R} \right\rbrack} - {\frac{1}{2}\overset{\_}{\beta}{\overset{\_}{Z}}^{\prime}\overset{\_}{\Lambda}\quad\overset{\_}{Z}}}$$H_{\alpha} = {\sum\limits_{n = 1}^{N}{\gamma_{{\alpha 0}\quad n}\left\lbrack {{\log\left( {{\overset{\_}{\alpha}}_{n}/{\overset{\_}{\alpha}}_{0n}} \right)} - \left( {{\overset{\_}{\alpha}}_{n}/{\overset{\_}{\alpha}}_{0n}} \right) + 1} \right\rbrack}}$$H_{\lambda} = {\sum\limits_{n = 1}^{N}{\gamma_{{\lambda 0}\quad n}\left\lbrack {{\log\left( {{\overset{\_}{\lambda}}_{n}/{\overset{\_}{\lambda}}_{0n}} \right)} - \left( {{\overset{\_}{\lambda}}_{n}/{\overset{\_}{\lambda}}_{0n}} \right) + 1} \right\rbrack}}$$H_{\tau} = {\gamma_{\tau 0}\left\lbrack {{\log\left( {\overset{\_}{\tau}/{\overset{\_}{\tau}}_{0}} \right)} - \left( {\overset{\_}{\tau}/{\overset{\_}{\tau}}_{0}} \right) + 1} \right\rbrack}$

In the second embodiment also, the posterior distribution P(J|B)corresponds to Q_(J) that maximizes F(Q), and the marginallog-likelihood log(P(B)) corresponds to the maximized F(Q).

Through the above described procedure, it is possible to calculate themarginal log-likelihood log(P(B)) of a virtual curved surface of acertain depth. Other procedures are the same as those of the firstembodiment, and therefore, description thereof will not be repeated.

[Meta Parameters of Hierarchical Prior Distribution]

Selection of meta parameters in the calculation of hierarchical priordistribution will be described.

When observation data obtained by other observation method is notavailable, the values of meta parameters are selected commonly,regardless of the lattice points, as in the first embodiment.

Specifically, bar τ₀, Γ_(β0) and γ_(τ0) can be selected in the samemanner as in the first embodiment. Further, meta parameters bar α_(0n),γ_(α0n), bar λ_(0n) and γ_(λ0n) can be represented as{overscore (α)}_(0n)={overscore (λ)}_(0n)={overscore (α)}₀γ_(α0n)=γ_(λ) _(0n)=γ_(α0)

Here, bar α₀ and γ_(α0) can be selected in the same manner as in thefirst embodiment.

(Third Embodiment)

In the first and second embodiments, the methods of estimating braincurrent sources have been described, assuming that observation dataobtained by other observation method is not available.

Recently, however, various methods have been developed to observe neuralactivities in the brain, and it has been made possible to obtainobservation data using a plurality of observing means under the sameexperimental conditions. Therefore, as the third embodiment,modifications from the first and second embodiments will be described,where observation data obtained by other observation methods such as MI,fMRI or PET are available.

First, when MRI observation data for the same person is available, it ispossible to obtain information related to the position of cerebralcortex or positions of other portions of the brain from MRI images. Insuch a case, it is possible to have the shape of virtual curved surfaceassumed in the first and second embodiments conforming to the shape ofthe cerebral cortex or the shape of other portions of the brain.

Further, when it is confirmed that a neural activity under certainexperimental conditions takes place in the cerebral vertex, it ispossible to omit searching in the depth direction and to estimatecurrent distribution by placing lattice points only on the cerebralcortex.

By fMRI or PET, it is possible to obtain information related to thelocation of neural activities. One must be careful, however, in handlingsuch information, because what is measured by fMRI is the amount ofblood flow and what is measured by PET is the amount of radio isotope.Though these are considered to increase at portions having higher neuralactivities, what is measured is not the current intensity derived fromthe neural activities.

It is considered that the degree of activities measured by f and PET areco-related to some extent to the current intensity derived from neuralactivities. It is not clear, however, whether the correlation is linearor non-linear. Further, MEG and EEG have temporal resolution in theorder of milliseconds, while fMRI has temporal resolution in the orderof seconds, and PET in the order of several tens of seconds.Specifically, what is obtained by fMRI or PET is an average over a longperiod of time of the measured results of MEG or EEG.

Dale et al. proposes a method of utilizing the result of observation ofMRI, fMRI or PET for the current source estimation by MEG or EEG inReference: A. M. Dale et al., Neuron. 26, pp.55-67 (2000). In thismethod, however, the activity information of fMRI or the like isdirectly applied as variance information of current.

From the view point of Bayesian estimation in accordance with thepresent invention, this corresponds to direct designation of the valueof inverse variance hyper parameter α_(n). As described above, however,the information obtained by f or the like is not related directly to thecurrent, and from the temporal view point, it is an average over a longperiod of time of the activities measured by MEG or EEG.

In contrast, in accordance with the present invention, the informationobtained by fMRI or the like is not directly designated as the varianceinformation of current but designated at the level of meta parameters ofhierarchical prior distribution. Meta parameters γ_(α0n) and γ_(λ0n) aremeta parameters that control reliability of information given by thehierarchical prior distribution. Specifically, the smaller the valuesγ_(α0n) and γ_(λ0n), the lower the reliability, and the higher thevalues γ_(α0n) and γ_(λ0n), the higher the reliability.

As described with reference to the second embodiment, the hierarchicalprior distribution in which localized condition and continuous conditionare combined is the same as assuming β⁻¹(A⁻¹+W·Λ⁻¹·W′) as the covariancematrix of current of the prior distribution, where A and Λ are diagonalmatrixes having {α_(n)|n=1, . . . , N} and {λ_(n)|n=1, . . . , N} asdiagonal components, respectively.

Specifically, when the value α_(n) is smaller, the current comes to havelarger variance, and when the value λ_(n) becomes smaller, thecovariance component derived from the continuous condition of currentbecomes larger. Meta parameters bar α_(0n) and bar λ_(0n) representexpected values of α_(n) and λ_(n) when there is no observation dataobtained from MEG (or EEG).

When it is expected that neural activity is vigorous at a currentlattice point n from fMRI or other observation method, the values ofmeta parameters bar α_(0n) and bar λ_(0n) are made smaller, and biasinformation may be entered to facilitate the current at this currentlattice point to assume a larger value.

In this manner, by introducing measured amount or information obtainedby other observation method having different temporal scale at the levelof meta parameters of hierarchical prior distribution, it becomespossible to enter information of other observation means as ambiguousinformation.

(Simulation Results)

Simulation results in accordance with the second and third embodimentswill be described.

In the following, an area including a certain visual cortex at a backportion of cerebral cortex obtained by MRI image of a person is used asa virtual curved surface, and the search in the depth direction isomitted.

As for the simulation data, currents having the current intensity of(1.0, 1.0, 0.5) (arbitrary unit) and currents having the currentintensity of (1.0, 0.0, 0.5) are caused to flow to areas V1, V2 and V5in the visual cortex, magnetic fields formed by these currents on theMEG sensor are calculated, and Gaussian noise having the S/N ratio of0.1 is added thereto, to provide the observation data.

FIGS. 28 and 29 represent results of Bayesian estimation consideringlocalized condition only. In the figures, a vertex of each trianglerepresents a lattice point prepared on the virtual curved surface, and acurrent dipole is allocated to each lattice point. The interval betweeneach of the lattice points is about 3 mm in average, and spatialresolution is high. Further, the virtual curved surface has acomplicated three dimensional shape as it is a separated part ofcerebral cortex. The figures are two-dimensional projection of the threedimensional shape. In the figures, reference characters V1, V2 and V5denote the areas V1, V2 and V5 of the visual cortex. In the figures,black circles represent points having high estimated current intensity.

FIGS. 30 and 31 represent result of Bayesian estimation considering bothcontinuous condition and localized condition. The results of estimationrecover the position and expansion of current distribution of almosttrue values. As can be seen from FIGS. 28 and 29, when Bayesianestimation is made using localized condition, the position of thecurrent is correctly estimated, whereas the expansion of the current isestimated to be too small. From this simulation, it is appreciated thatintroduction of continuous condition is effective when spatialresolution is high.

Next, a simulation was made with the S/N ratio of noise increased to 0.2and the current distribution unchanged. In that case, even whenhierarchical prior distribution combining the localized condition andcontinuous condition was used, the expansion of current distributioncould not be recovered, because of the influence of severe noise.

Thereafter, a simulation was made assuming that information ofobservation data obtained by other observation means is available underthe bad condition as above. In the simulation, it was assumed thatinformation indicating that activities in areas V1, V2 and V5 arevigorous is obtained from the observation data of fMRI. Considering thefact that the active areas of fMRI do not always correspond to currentactivities and that fMRI corresponds to an average over a long period oftime of MEG observed data, it was assumed that the areas indicating highactivities in fMRI were areas V1, V2 and V5 expanded by about 20 mm inradius. The same fMRI information was used when the current intensitiesat V1, V2 and V5 were (1.0, 1.0, 0.5) and (1.0, 0.0, 0.5), respectively.

By the setting above, a situation is simulated in which the active areaof fMRI does not always correspond to the current intensity, though itis correlated to some extent. Under such setting, meta parameter oflattice points corresponding to the active area of fMRI was set to barα_(0n)=bar λ_(0n)=2×Tr(G′·G)/tilde N, and meta parameter of otherlattice points was set to bar α_(0n)=bar λ_(0n)=50000×Tr(G′·G)/tilde N.Here, the term “tilde” preceding a variable means that the variable hasthe sign “˜” thereabove.

The values of γ_(α0n) and γ_(λ0n) were commonly set toγ_(α0n)=γ_(λ0n)=0.1 for all points.

FIGS. 32 and 33 represent the results of estimation at this time. Fromthese figures, it can be understood that even when there is considerablenoise, correct estimation is possible by adding information obtained byother observation means.

Although the present invention has been described in detail, it isclearly understood that the same is by way of illustration only and isnot to be taken by way of limitation, the spirit and scope of thepresent invention being limited only by the terms of the appendedclaims.

1. A brain current source estimating method for estimating, based on anelectromagnetic field observed outside a scalp, a position of a currentsource as a source of said electromagnetic wave existing in the brain,comprising the steps of setting, in the brain, a plurality of virtualcurved surfaces having depths from brain surface different from eachother and shapes not intersecting with each other, and setting latticepoints on each of said virtual curved surfaces; estimating, on each ofsaid virtual curved surfaces, a: current distribution for recoveringsaid observed electromagnetic field; and based on an expansion of thecurrent distribution estimated on said virtual curved surfaces and adifference between the electromagnetic field recovered based on saidcurrent distribution and said observed electromagnetic field,identifying a virtual curved surface at which said expansion and saiddifference attain relative minimums among said plurality of virtualcurved surfaces as a true curved surface on which said current sourceexists.
 2. The brain current source estimating method according to claim1, wherein said step of estimating said current distribution includesthe step of determining posterior probability by Bayesian estimationmethod from prior distribution and observation data of saidelectromagnetic field; and said step of identifying as a true curvedsurface on which said current source exists includes the step ofidentifying a virtual curved surface of which corresponding saidposterior probability attains the maximum, among said virtual curvedsurfaces.
 3. The brain current source estimating method according toclaim 2, wherein said step of estimating a current distribution includesthe step of identifying a first virtual curved surface closest to saidbrain surface and having posterior probability attaining a relativemaximum, among said plurality of virtual surfaces, while successivelymoving from a virtual curved surface on the side of the brain surface toa deeper side; and said step of identifying a curved surface as a truecurved surface on which said current source exists includes the steps ofidentifying a localized current distribution corresponding to a point ofrelative maximum of said current distribution, on said first virtualcurved surface, separating a plurality of local surfaces each includingsaid localized current distribution, and fixing, among said plurality oflocal surfaces, local surfaces other than a local surface as an objectof identification, moving said local surface as an object ofidentification in the depth direction, and identifying positions wheresaid posterior probability attains the relative maximum, successivelyfrom the side closer to said brain surface.
 4. The brain current sourceestimating method according to claim 3, wherein in said step ofestimating a current distribution, said current distribution isestimated with a first spatial resolution; said method furthercomprising the step of re-estimating said current distribution with asecond spatial resolution higher than said first resolution andresolution of said plurality of virtual curved surfaces in the depthdirection being improved.
 5. The brain current source estimating methodaccording to claim 1, wherein said step of estimating a currentdistribution includes the step of setting, when said currentdistribution is estimated in accordance with Bayesian estimation, ahierarchical prior distribution in said Bayesian estimation usingobservation data obtained by other observation method independent ofsaid observation of electromagnetic field for said estimation of thecurrent source.
 6. A program for a computer for estimating, based on anelectromagnetic field observed outside a scalp, a position of a currentsource as a source of said electromagnetic wave existing in the brain,to have the computer execute the steps of: setting, in the brain, aplurality of virtual curved surfaces having depths from brain surfacedifferent from each other and shapes not intersecting with each other,and setting lattice points on each of said virtual curved surfaces;estimating, on each of said virtual curved surfaces, a currentdistribution for recovering said observed electromagnetic field; andbased on an expansion of the current distribution estimated on saidvirtual curved surfaces and a difference between the electromagneticfield recovered based on said current distribution and said observedelectromagnetic field, identifying a virtual curved surface at whichsaid expansion and said difference attain relative minimums among saidplurality of virtual curved surfaces as a true curved surface on whichsaid current source exists.
 7. The program according to claim 6, whereinsaid step of estimating said current distribution includes the step ofdetermining posterior probability by Bayesian estimation method fromprior distribution and observation data of said electromagnetic field;and said step of identifying as a true curved surface on which saidcurrent source exists includes the step of identifying a virtual curvedsurface of which corresponding said posterior probability attains themaximum, among said virtual curved surfaces.
 8. The program according toclaim 7, wherein said step of estimating a current distribution includesthe step of identifying a first virtual curved surface closest to saidbrain surface and having posterior probability attaining a relativemaximum, among said plurality of virtual surfaces, while successivelymoving from a virtual curved surface on the side of the brain surface toa deeper side; and said step of identifying a curved surface as a truecurved surface on which said current source exists includes the steps ofidentifying a localized current distribution corresponding to a point ofrelative maximum of said current distribution, on said first virtualcurved surface, separating a plurality of local surfaces each includingsaid localized current distribution, and fixing, among said plurality oflocal surfaces, local surfaces other than a local surface as an objectof identification, moving said local surface as an object ofidentification in the depth direction, and identifying positions wheresaid posterior probability attains the relative maximum, successivelyfrom the side closer to said brain surface.
 9. The program according toclaim 8, wherein in said step of estimating a current distribution, saidcurrent distribution is estimated with a first spatial resolution; saidmethod further comprising the step of re-estimating said currentdistribution with a second spatial resolution higher than said firstresolution and resolution of said plurality of virtual curved surfacesin the depth direction being improved.
 10. The program according toclaim 6, wherein said step of estimating a current distribution includesthe step of setting, when said current distribution is estimated inaccordance with Bayesian estimation, a hierarchical prior distributionin said Bayesian estimation using observation data obtained by otherobservation method independent of said observation of electromagneticfield for said estimation of the current source.
 11. A brain currentsource estimating apparatus for estimating, based on an electromagneticfield observed outside a scalp, a position of a current source as asource of said electromagnetic wave existing in the brain, comprising:virtual curved surface setting means for setting, in the brain, aplurality of virtual curved surfaces having depths from brain surfacedifferent from each other and shapes not intersecting with each other,and setting lattice points on each of said virtual curved surfaces;current distribution estimating means for estimating, on each of saidvirtual curved surfaces, a current distribution for recovering saidobserved electromagnetic field; and current source identifying means foridentifying, based on an expansion of the current distribution estimatedon said virtual curved surfaces and a difference between theelectromagnetic field recovered based on said current distribution andsaid observed electromagnetic field, a virtual curved surface at whichsaid expansion and said difference attain relative minimums among saidplurality of virtual curved surfaces as a true curved surface on whichsaid current source exists.
 12. The brain current source estimatingapparatus according to claim 11, wherein said current distributionestimating means includes posterior probability determining means fordetermining posterior probability by Bayesian estimation method fromprior distribution and observation data of said electromagnetic field;and said current source identifying means includes virtual curvedsurface identifying means for identifying a virtual curved surface ofwhich corresponding said posterior probability attains the maximum,among said virtual curved surfaces.
 13. The brain current sourceestimating apparatus according to claim 12, wherein said currentdistribution estimating means includes shallowest virtual curved surfaceidentifying means for identifying a first virtual curved surface closestto said brain surface and having posterior probability attaining arelative maximum, among said plurality of virtual surfaces, whilesuccessively moving from a virtual curved surface on the side of thebrain surface to a deeper side; and said current source identifyingmeans includes localized current distribution identifying means foridentifying a localized current distribution corresponding to a point ofrelative maximum of said current distribution, on said first virtualcurved surface, local surface extracting means for separating aplurality of local surfaces each including said localized currentdistribution, and local surface position identifying means for fixing,among said plurality of local surfaces, local surfaces other than alocal surface serving as an object of identification, moving said localsurface as an object of identification in the depth direction, andidentifying positions where said posterior probability attains therelative maximum, successively from the side closer to said brainsurface.
 14. The brain current source estimating apparatus according toclaim 13, wherein said current distribution estimating means estimatessaid current distribution with a first spatial resolution and thereafterre-estimates said current distribution with a second spatial resolutionhigher than said first resolution and resolution of said plurality ofvirtual curved surfaces in the depth direction being improved.
 15. Thebrain current source estimating apparatus according to claim 11, whereinsaid current distribution estimating means includes means for setting,when said current distribution is estimated in accordance with Bayesianestimation, a hierarchical prior distribution in said Bayesianestimation using observation data obtained by other observation methodindependent of said observation of electromagnetic field for saidestimation of the current source.
 16. A recording medium recording acomputer program for estimating, based on an electromagnetic fieldobserved outside a scalp, a position of a current source as a source ofsaid electromagnetic wave existing in the brain, said program comprisingthe steps of: setting, in the brain, a plurality of virtual curvedsurfaces having depths from brain surface different from each other andshapes not intersecting with each other, and setting lattice points oneach of said virtual curved surfaces; estimating, on each of saidvirtual curved surfaces, a current distribution for recovering saidobserved electromagnetic field; and based on an expansion of the currentdistribution estimated on said virtual curved surfaces and a differencebetween the electromagnetic field recovered based on said currentdistribution and said observed electromagnetic field, identifying avirtual curved surface at which said expansion and said differenceattain relative minimums among said plurality of virtual curved surfacesas a true curved surface on which said current source exists.
 17. Therecording medium according to claim 16, wherein said step of estimatingsaid current distribution includes the step of determining posteriorprobability by Bayesian estimation method from prior distribution andobservation data of said electromagnetic field; and said step ofidentifying as a true curved surface on which said current source existsincludes the step of identifying a virtual curved surface of whichcorresponding said posterior probability attains the maximum, among saidvirtual curved surfaces.
 18. The recording medium according to claim 17,wherein said step of estimating a current distribution includes the stepof identifying a first virtual curved surface closest to said brainsurface and having posterior probability attaining a relative maximum,among said plurality of virtual surfaces, while successively moving froma virtual curved surface on the side of the brain surface to a deeperside; and said step of identifying a curved surface as a true curvedsurface on which said current source exists includes the steps ofidentifying a localized current distribution corresponding to a point ofrelative maximum of said current distribution, on said first virtualcurved surface, separating a plurality of local surfaces each includingsaid localized current distribution, and fixing, among said plurality oflocal surfaces, local surfaces other than a local surface serving as anobject of identification, moving said local surface as an object ofidentification in the depth direction, and identifying positions wheresaid posterior probability attains the relative maximum, successivelyfrom the side closer to said brain surface.
 19. The recording mediumaccording to claim 18, wherein in said step of estimating a currentdistribution, said current distribution is estimated with a firstspatial resolution; said method further comprising the step ofre-estimating said current distribution with a second spatial resolutionhigher than said first resolution and resolution of said plurality ofvirtual curved surfaces in the depth direction being improved.
 20. Therecording medium according to claim 16, wherein said step of estimatinga current distribution includes the step of setting, when said currentdistribution is estimated in accordance with Bayesian estimation, ahierarchical prior distribution in said Bayesian estimation usingobservation data obtained by other observation method independent ofsaid observation of electromagnetic field for said estimation of thecurrent source.